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ABSTRACT 



A computer model for axially symmetric dielectric radomes in the near field of 
a circular aperture has been developed. This model, based on a method of 
moments solution for bodies of revolution, is used to determine the electric field 
pattern of a linearly polarized circular aperture radiating through a dielectric 
radome. The program is written in FORTRAN and can accommodate rotationally 
symmetric radome shapes. Arbitrary illumination distributions for the aperture can 
be specified, and phased scanning of the aperture in both azimuth and elevation 
is allowed. Computational run time has been reduced by taking advantage of 
mode symmetry. A separate program has been developed to compute gain, and 
together these programs are used to determine the pattern degradation and gain 
loss due to the presence of a radome in the near field. 
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I. INTRODUCTION 



A radome is a structure designed to protect an antenna in its operating 
environment. Aircraft radomes must protect radar antennas from high 
aerodynamic forces while minimizing drag, weight, and interference on its 
electrical operation. The presence of a radome can affect the gain, beamwidth, 
sidelobe level, and direction of boresight, as well as change the VSWR and 
antenna noise temperature. This thesis involves the continued development of a 
computer code that predicts antenna radiation patterns looking through 
aerodynamic radomes. 

The ability to accurately predict the effect of a radome on the operation of an 
antenna early in the design process is desirable. Typically, designers of the 
earliest radomes were concerned more with structural considerations than 
electrical performance. However, after the near failure of the U.S. Army's 
"Pathfinder" program of radar-equipped B-24 "Liberator" bombers due to poor 
electrical performance of their radomes, it became clear that," There was a serious 
need for developing sound theoretical procedure for the electrical design of 
radomes." [Ref. 1] For most practical radomes the surface geometry is too 
complex to yield closed-form mathematical expressions describing their 
performance. It is possible to predict radome performance using numerical 
solutions of Maxwell's equations. The program described in this thesis is based 
on a numerical technique called the method of moments (MM). The particular 
form of the solution used here is based on a formulation by Mautz and Harrington 
[Ref. 7] for bodies of revolution. Francis [Ref. 2] developed a computer code 
LDBORMM.F that specialized the Mautz and Harrington solution to ogive- 
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shaped radomes. There are no restrictions on the location of the radome with 
respect to the antenna, and therefore near-field radomes can be analyzed. 

The goal of this thesis was to extend the capabilities of this program. 
Significant improvements include: taking advantage of symmetry to reduce 
computation time, applying realistic radar antenna patterns for comparison to 
experimental results, and the computation of gain loss due to the radome. 

A. DESCRIPTION OF THE PROBLEM 

The basic system analyzed is a radome consisting of a body of revolution 
(BOR) located coaxial with a circular antenna, as in Figure 1.1. The BOR 
constraint is satisfied by many aerodynamic radomes and greatly simplifies the 
MM solution , which will be described in more detail in following chapters. The 




Axis of 
Symmetry 



ogive shape of the radome in Figure 1.1 is the primary one of interest, but the 
program can also model other bodies of revolution including: a cone, disk, 
hemisphere, and parabola, as shown in Figure 1.2. These were primarily added for 
validation purposes. 
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There are no constraints on the size, shape and location of these objects 
provided they are coaxial with the antenna. However, because MM requires the 




solution of a matrix equation, computer run time and memory will limit the size of 
the bodies that can be analyzed. A unique feature of this program is that it takes 
into consideration the fact that the radome may be in the near field of the antenna 
as in Figure 1.3. The antenna aperture used in the program is circular, but the 
program can be easily altered to accommodate other shapes. There are no 
constraints on the aperture size, scan angle, or polarization. 

The remainder of this thesis is divided into three major sections. Chapter II 
covers theoretical development of the problem and method of solution. The MM 
solution is described in Chapter II along with a discussion of mode symmetry and 
how it can be used to advantage in obtaining the radome current. Chapter II also 
describes how the total field is calculated from the currents. Chapter III covers 
the programs used for generating patterns based on the methods of 
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Chapter II. The first part describes the program RADOME.F which calculates the 
current and electric field. The second part describes the program GAIN.F which 




Figure 1.3 Radome Located in Near Field of Antenna 

calculates the gain of the system. The final chapter discusses the results and 
presents several conclusions based on the calculated radome patterns. 
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II. ANALYSIS 



A. GENERAL DESCRIPTION OF ANALYSIS 

The objective of the analysis is to determine the radiation pattern of the 
antenna and its gain in the presence of a radome. In order to analyze the system, 
the spherical coordinates of 1 lgure 2.1 are used. The total electric field in the 



total electric field is the sum of the electric field due to the currents on the 
antenna (radiated field), and the electric field due to the currents on the radome 
(scattered field). The procedure is to assume a current distribution on the 
antenna, then use the MM to determine currents induced on the radome. In 
practice, scattering from the radome back toward the antenna will cause a change 
in the initial current distribution. This effect will be small for a well-designed 
radome and is ignored in the analysis. Scattering back toward the antenna will 
primarily affect the VSWR looking into its terminals. 

1. Determination of Gain 

Gain is a measure of the ability of an antenna to concentrate radiated 
power in a particular direction. The maximum value of gain for a lossless antenna 
is the directivity Go, which is the ratio of the maximum radiation intensity to the 
average radiation intensity [Ref. 3:p 610], 



direction of unit vectors |0.p) is required to compute the gain in the far field. The 




§u(e,n>)d a 



( 2 - 1 ) 



s 
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where U[d,o) is the radiation intensity and dQ. = Sin ddOdv is the differential 
solid angle. Throughout this paper a lossless antenna will be assumed and the 




Figure 2.1 Coordinate System 



terms gain and directivity used interchangeably. 

Radiation intensity is the time-averaged power per unit solid angle, and 
is proportional to the magnitude of the electric field squared 



U = 




( 2 - 2 ) 



where 
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(2-3) 



E 



E a 




and E J is the electric field from the antenna, E s is the electric field scattered by 
the currents induced on the radome. and r| is the impedance of free space. 

2. Determining the Electric Field 

The total electric field is calculated from the current distributions on the 
antenna and radome as illustrated in Figure 2.2. A current distribution J a on the 
antenna is assumed known. The currents induced on the radome are dependent 
on the current distribution on the antenna as shown in Figure 2.3, and this is 
taken into account in the method of moments solution. 

Maxwell's equations, are used to determine the electric field at a point in 
space due to a current density J. When both the radome and antenna are 
present, the total current density is J = J a + J s . A solution of Maxwell's 
equations is 

E = -jwR Lv(v.r), (2-4) 

cope ' ’ 

where 

n = -^-ffj(r.e,(J))— ds', (2-5) 

4k JJ r 

and 



r = 



R-R' 



( 2 - 6 ) 



In (2-4) R is the vector potential, R the position vector to the observation point 
and R' the position vector to a source point on the antenna. 
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Figure 2.2 Electric Field from Currents on Antenna and Radome 




Figure 2.3 Electric Field Impinging on Radome 
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B. ANALYTICAL TECHNIQUE 

Closed form mathematical solutions of the gain and radiation patterns are only 
possible for a limited number of applications that neglect the presence of the 
radome. In general, determining the radiation pattern and gain with the radome 
requires a numerical solution of Maxwell’s equations. The radome geometry 
results in an unknown current distribution so (2-4) cannot be solved directly. 
However, boundary condition ^ can be applied on the radome surface which yield 
an integral equation for the current. The integral equation will be solved 
numerically using the method of moments. 

1. E-field Integral Equation and the Method of Moments 

The method of moments (MM) is a procedure used to solve an integral 
equation obtained by enforcing (2-4) on the radome surface [Ref. 5], Equation (2- 
4) is used to determine the electric field E a incident upon the surface of the 
radome from the currents on the antenna J a . Initially the surface will be assumed 
a perfect electric conductor (PEC), in which case the total E tangent to the surface 
of the radome is equal to zero. Later a correction term will be added to the perfect 
conductor solution to account for the dielectric properties of the radome. 



The boundary condition on the PEC radome is 




(2-7) 



As applied to (2-4), 




( 2 - 8 ) 



where 
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(2-9) 



a~Jkr 

G = - , 

4^r 



and 



r = 



R-R'. 



( 2 - 6 ) 



Equation (2-8) is the E-field integral equation (EFIE). The only unknown is J s , 
and once determined it is possible to calculate the scattered field by applying (2-4) 
to the radome. 

To solve (2-8), the current J s is expanded into a series, 

j, = a- »o) 

1=1 

The Jj are basis functions, and the I, are the expansion coefficients to be 
determined, 

Mautz and Harrington applied this technique to the specific case of 
surfaces S generated by revolving a plane curve about the z-axis [Ref. 5], The 
plane curve is approximated by a series of straight-line segments connecting N p 
points. The surface has a circular cross section in the azimuth variable thus the 
body is represented by a series of "faceted rings" as shown in Figure 2.4. 
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The surface coordinate system is cylindrical [p,o,Z 1, and t is a vector 

tangent to the surface and orthogonal to o at every point. In order for the series 
expansion of the current (2-10) to approximate an arbitrary surface current density 
J s , independent sets of functions are defined as follows: 




Figure 2.4 Segmented Radome with Current Basis Functions 
Depicted 



in t: J' nj 


- Tj it) 

- t J ' e jno 

Ps 


n = 0, ±1, 


• • > ±°° 


y = 1,2, 


..,Np -2 


in 0\ J a : 

nj 


= 0 w e jno 

Ps 


n = 0,±1, 


• .±°° 


j= 1,2,.. 


£ 

1 

i 



Tj(t) is a triangle function extending over two segments, and Pj[t j is a 

pulse function extending over a single segment. For an explanation ot these 
functions see [Ref. l:pp 16-18]. Each value of n is referred to as an azimuthal 
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mode. The series expansion of (2-10) in terms of the chosen basis functions 
becomes. 



•*.= I 



sp-2 _ Np-I 

I '1A + 1 '„’A 

L 7=1 7=1 



( 2 - 12 ) 



In order for the series to be an exact solution, the sum must be extended over 
infinite modes -°o</7<oo. Of course the series must be truncated, and the 
maximum number of modes for convergence n mait is dependent upon the antenna 
scan angle and radome size. 

When the series expansion is applied to the integral equation (2-8), 



j Sp- 2 

= I It SI 



n=-°° 7=1 S, 
+~ \p-l 






E- = I I'Jj 

" = — 7=1 S, . 



jco^ 6 nj G - 



coe 

zl 

coe 



(v-.j^jvc 



ds' 

ds'. 



(2-13) 



In order to solve for the unknowns l' nj and l ° nj , both sides of (2-13) are multiplied 

by testing functions and UJ ° mi , and integrated. Using Galerkin's Method 

[Ref. 8], the testing functions are chosen to be the complex conjugates of the basis 
functions J /,e , 




g~J m 0 



g-jm* 



(2-14) 
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Multiplying both sides of (2-13) by each of the testing functions and 
integrating results in a set of 2 N p - 3 simultaneous linear equations for each value 



of n In matrix form, these equations are written as 



U 




( 2 - 15 ) 



where 




is the excitation vector. 




is an impedance matrix and 



contains the 



unknown current coefficients. Submatrices associated with the t. and 0 
components (2-15) can be identified as follows: 





u r 








\z" 1 

n 




r° 

n 






[I'n 






U 9 

n 








Z°' 

n 




£00 

n 






I 9 

n 





n - 0,±1,±2,. 



1 * 



( 2 - 16 ) 



The unknown coefficients are determined by solving the matrix equation 



I 



-i 



U 



( 2 - 17 ) 



The MM solution for the unknown current coefficients l' nj and l° nj in vector form 
is 













\z" 1 

*- n 




[ Z "l 




-1 










r 








Z ol 




£00 








u 9 






n 
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n 




n 






_ 


n 


_ 



( 2 - 18 ) 



Once these are determined, it is possible to compute the electric field due to the 

current density distribution J s on the surface of the radome at any point in space. 
The elements of each submatrix of the excitation vector II are 
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( 2 - 19 ) 



U‘ = — ff UJ 1 • E J 'ds 

nt JJ w ni "-tan 

s 

U° =±\[LU° • E a f/s, 

m „ JJ ni <an ’ 



and the elements of each of the submatrices of the impedance matrix 




are 






yto 



yot 

n 



yoo 

n 



s, 5, L 

S/ 5/ L 

s / */ L 



ds'ds 
ds'ds 
ds'ds 
ds'ds . 



( 2 - 20 ) 



The fact that the testing functions U ) ^ are orthogonal in q to due to 
the presence of the exponentials has been used to eliminate elements of (2-20) 
with n * m [Ref. 6]. This allows each mode n to be treated independently, 
thereby simplifying the matrix inversion problem. Thus, by constraining the 
surface to be rotationally symmetric, the azimuthal dependence of the current is 
expressed as a sum over exponentials (a Fourier series) and each mode n is 
independent. 



a. Impedance of Thin-Shell Dielectric Radomes 



The impedance matrix 



of the MM solution assumes that the 
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radome surface is a PEC, where (2-7) is applied as a boundary condition of (2-4). 
For thin-shell dielectric radomes, a correction term is added to the impedance 
matrix that allows the modeling of radomes of arbitrary dielectric constant [Ref.9]. 

When the electromagnetic wave from the antenna is incident on a 
dielectric radome that has an intrinsic impedance 



n = -l^ 



( 2 - 21 ) 



different from that of free space (b 0 ). there will be reflection and refraction as 
shown in Figure 2.5 [Ref. 3:p 397], The total fields are 



E = V + V 

H =H S +H a 



( 2 - 22 ) 



These fields satisfy Maxwell's equations 



V x ( E a + E* ) = -jo)fd 0 (fi a + H s ) 



( 2 - 23 ) 



and. 



V x (H a + H s ) = jcoe 0 [V + E s ) + jco[e - e 0 )(E a + 



( 2 - 24 ) 



where the second term exists only when there is material present with a 
permittivity different than free space. This term is a polarization current J. where 



J S j0)(e- £ 0 )E. 



( 2 - 25 ) 
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Maxwell's equations are also satisfied by the field from the antenna |E a ,H a j. 
Since I^H) and ( E a , H a ) satisfy Maxwell's equations, by supeiposition and (2-22) 
(E S ,H S ) also satisfy Maxwell's equations. Therefore 



VxH s = jcoV + J 



(2-26) 




Figure 2.5 Discontinuity of Permittivity for Dielectric 
Radomes 



For a radome consisting of a thin dielectric shell, the polarization 
current has a large tangential component and small normal component. If only the 
tangential component is significant, then (2-8) can be modified to include 
dielectric properties of the radome as follows: 



16 



( 2 - 27 ) 



“tan 



ja)(e ~ 



r r “* / 


r i 


ff J s Gds + ^-V 


V • 1 1 J ,Gds 


j coe 


JJ 

s 



The first term on the right Nide is the impedance loading term; the integro- 

differential portion remains unchanged from the PEC assumption. 

After multiplying both sides of (2-27) by the testing functions 
(liC, HC) and integrating, the impedance matrix |zj of (2-15) becomes 



[Z] = [Z MM ] + [Z L ], 



( 2 - 28 ) 



where the elements of 



-MM 



are given by (2-20), and 



Z "»L J =S\‘H!'* J nZ L ds 



-/to 
**L n 



= jf/^.J v „Z,ds = 0 , (by orthogonality of t, and o) 

S 

= |] ui; • j '„ n z L ds = a , (bv orthogonality of o. and t) 

S 

Z“] 



-jot 

*-Ln 



( 2 - 29 ) 



and 



Zl = 



1 

Ja)(e ~ f 0 )f 



( 2 - 30 ) 



The thickness t in the impedance loading term is necessary to convert the volume 
current density J, of Figure 2.5 to a surface current density J s . 
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Solving (2-17) yields 



I 





(2-31) 



The radome impedance Z L can be rewritten in terms of dielectric constant e r , 
thickness to wavelength ratio n . and the loss tangent tan 8 [Ref. 2:pp 45-46], 





60 




j(e r -1) + e r tan 5 



(2-32) 



Note that both (2-29) and (2-32) allow the modeling of radomes whose 
permittivity may be complex, e = e' - je" . 
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2. Mode Symmetry 

When the surface current density J s is expanded into a series as in (2- 
12), the sum is extended over azimuthal modes n = 0,±l,±2, ,±n maK . The 
mode index appears in the exponential term of basis functions (2-1 1) and 
testing functions ZZ/Jv (2-14). The exponential factor can be written as, 

e jn = cos{no) + jsinino). (2-33) 



For each n there is a current vector l„ given by (2-18). Mode symmetry refers to 
a relationship between the positive and negative mode vectors l„ and l_„ . If mode 
symmetry exists, then only one matrix equation needs to be evaluated rather than 
two. This greatly reduces the computational effort. 



Recall the definition of the excitation vector 



0 



u'„, 

5 

U =±\\LU a *t°dS (2-19) 

m n } J ni ten v 



The electric field from the antenna is 

V =E a „(R, e,o)R + E a e (R,d, <t>)e + £*(ff, 0 , 0 ) 0 , (2-34) 

which on the radome surface can be expressed in terms of the two orthogonal 
components t and 0 



e* =£■;(!,*)» +c; ( 1 , 0)1 



(2-35) 
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In the far field of the antenna, the R component in (2-34) decays to zero and Ef 
consists only of E*. Therefore as shown in Figure 2.6, 



E?(t,o) = Ef(t)cos0 

E*(t,o) = E°{t)sin<t>. (2-36) 



x 




Figure 2.6 Components of the Antenna Electric-Field on the Radome 



Plugging these far-field components into (2-19) 



U' ni = J j coso cos(/7tf>) + jsin(np)]d0 

J = 0 Psl ' ' 



o=0 

2t 



U* m = j j Sino cos(/?<^) + Jsin(/7o)]t/0. (2-37) 

/=e ' 



0 = 0 
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Only the o integration is dependent on n Expanding the integral in o and 
explicitly including the sign o! n gives 

2 T 



j cos 

(5 = 0 



0 



c o s ( ± no ) + jsin 




2<t 



2t 



J cos o cos( ±no)f/p + J J cos p sinj ±no )do. (2-38) 



<> = 0 



0 = 0 



where 



2 - 

| cos os\n(±no)dd = 0. (2-39) 

0 = 0 



due to the fact that the integrand is an odd function of 0 , and the integration is 
over one period 0 < 0 < 2 n . By trigonometric identity 



2-t 2k 2 k 

J COS 0 COS( ±/70 jrfO = J COS (±/7-lJoko+ J COS (±/7 + l)p 

0=0 =0 o =0 



do. (2-40) 



The right hand side is only non zero if n = ±1. Furthermore, because cosine is an 
even function, the term U * ni is independent of the sign of n. Thus 

U[ nJ =+U[ nJ . (2-41) 

For the second equation of (2-^7). the sign of U° nj is dependent on the integral 
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(2-42) 



2k 

J sine cos(±/7o) + js\r\{±nomo = 



o = 0 
2k 



2 k 2 k 

J s\nocos[±no ] jdo + j J sin <z>sin(±/70 )tfo. 



0 = 0 



o = 0 



where 



2 - 

\ sinocos{±no\do = 0, 

o = Q 



(2-43) 



by symmetry. By trigonometric identity 



2 K 2 t 2t 

| sin 0 sin[±nojtfo = | cos \±n - ijopo - J cos(±n + l) 

o =0 0=0 0=0 



0 



do. (2-44) 



As can be seen in (2-44), the excitation vector U c nj is dependent on the sign of n 
due to the presence of the negative sign. Therefore 

U° =-U° . (2-45) 

-nj +nj ’ 



In the near-field of the antenna the 0 dependence is more complicated 

than the sine and cosine functions in (2-36). However, it has been verified 
numerically that E * and E a R are always even functions of o and E a 0 is always odd. 



This conclusion assumes that a symmetric amplitude and linear phase distribution 
excites the antenna. 



In light of the above results, the elements of the impedance matrix 



-MM 



of (2-18) have the following relationship between positive and negative modes: 
[Ref. 6] 
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[Z 



MM 



z' ' 

l MM (-n) 

z ot 



■ yt o 

-J9 0 

L MM f-n 



z' ' 

1 MM(+n 

-I” 

^MMf+n 



f-Z' 5 

ij [ t MM *n 

fz® 9 



ij 

'J 



(2-46) 



By inspection of (2-30), the elements of the load impedance matrix Z, are 
independent of the mode n. because /Z//„® and J'j n are complex conjugates of each 
other. Therefore, 
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(2-47) 



Rewriting (2-31) results in 



= [v][u 



[2-48) 



where 



Zmm + Z L 



i-i 



(2-49) 



By taking the inverse through partitioning [Ref. 8] it can be shown that the mode 
symmetry survives the inversion, and 
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12-50) 



By multiplying out (2-48), 
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Equations (2-41 ). (2-45), and (2-50) relate the positive and negative modes of (2- 
51), as follows: 



T = V" U 1 -V h [-U* ) = +1' 

•-n ~mj u ~nj +mj \ +nj ) + n 

V =-T' U[ n , +T° (-U* ,) = -\* . ( 2 - 52 ) 

~n - mj +mj + nj ] + n v ; 



This is an extremely important result. Using mode symmetry reduces the number 
of calculations required to find the current coefficients by nearly half in most 



cases. 



3. Measurement Matrices for the Scattered Field 

Once the current J s is obtained the scattered field E s can be determined 
by summing the contributions of the individual current elements in the far-field as 
shown in Figure 2.7. This is accomplished by use of a measurement vector with 
elements 



In (2-53), a unit radiated plane wave is weighted by the current element at each 
point on the surface. Neglecting the R component in the far-field, the 
measurement vector separated into 0 and P components is 




(2-53) 




In 

far-field 



Figure 2.7 Scattered Field 
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(2-54) 




e 



- jki hu +yu +luz 



( 0 . j S )ds\ 



where U . U . and UJ are the direction cosines. 



u = sinflcoso 
u = sin 0sin0 

uj - cos0. (2-55) 



Thus, the components of the scattered field at a far-field point are the 
superposition of all surface segments 



£S _ zjhll g-J kr 
*4 nr 



I Kl 




-Jk'i p-Jkr \ no / 

A nr C Zu n m 1 ™' 
m 



(2-56) 



4. Closed-Form Far-Field Antenna Pattern 

When computing the antenna contribution to the total field in the far 
zone the higher order j terms can be neglected, and the radiation integral can be 

reduced to a closed form. The source of this field is assumed to be a uniformly- 
excited circular-aperture antenna scanned to a direction (0 s ,0 s j. Uniformly 

excited refers to a constant amplitude and a linear phase distribution of current on 
the antenna. The /?, 6 and p components in the far-field are 



- 0 
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-jk e~ Jkr 

4 k r 
-jk e 
4/r r 



(L* + IN*) 

( l s + n^.). 



. where L 0 , /.* /V* and N e are given by 



(2-57) 



L = \\Me- Jkr cosv ds' = 0. (M = 0), 

S' 

N = JJ je' Jkr Q0S¥ ds' = j J (k J h + y + z J z y jkrC0S¥ ds' (2-58) 

S’ s' 



where r' and yt' are as shown in Figure 2.8 [Ref. 10:pp. 455], If the antenna is 
x-polarized, then N = N z = 0 and 






jj J H e JkrC0Sv ds'. 



(2-59) 



The equivalent in spherical coordinates is 

N„ = Jj J H cosOcosoe Jkr COi¥ ds' 

s 

N 0 = \\j H s\noe JkrcoS¥ 'ds'. (2-60) 

5 
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The phase of the source point relative to the field point, /"COS i//' , is a length that 
can be written in Cartesian coordinates as 



r'cos y/' = h u + y'u + z'uj, ( 2 - 61 ) 




Figure 2.8 Circular Aperture Antenna Source Point 



where u, u , and LU are the direction cosines (2-55). Including beam scanning 
adds an exponential phase factor to the current 



J 



H 



j e Jk\KU, + yv,*zw,) 



( 2 - 62 ) 
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where U s . U s and LU S are the direction cosines of the scan angles [d s ,o s \ which 

also obey (2-55). J 0 is the magnitude distribution, assumed to be constant, and 
[x',y', Z' | are the coordinates of a source point. For an antenna iying in the x-y 

plane, Z' = 0. Therefore, N 8 can be rewritten as 

k I i rr jk[H'u s +yu.) -jkix’u +u‘u\ , , 

N 8 = J 0 cos 6cos o J J e e ' y ds . 

s 

. rr jk[x'[u s -u\+u’[u t -u\\ , , ^ , 

= J 0 cos 6>cos ojj e ds . (2-63) 



By setting X' = p' cos o' and y' = p' Sin o', the exponent term of (2-63 ) become ; 



kp' coso'(Uj -u) + sin 0 '(z/ s -u 



(2-64) 



Defining the variables 



fl = (u,-u) 

B=(o, -o), (2-65) 

and considering these to be orthogonal vectors as in Figure 2.9, it is possible to 
rewrite (2-64) as 



kp'iff 2 + B 2 



R B 

coso— + sin 0 ' 



MR 2 + B 2 



Vfl 2 + B 2 



( 2 - 66 ) 
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From Figure 2.9 the factors containing R and B in the brackets can be replaced 
with cos / and sin/ 



kp'\R 2 + B 2 [coso'cosy + sinp'sin/ 



(2-67) 



By trigonometric identity. 



k, A R 2 + B 2 



COSI 



(®- - r) 



( 2 - 68 ) 




Figure 2.9 Representation of R and B as vectors 



Using this form of the exponent in the integral of (2-63), gives 



iV tf = J 0 COS8COSO 



J ]e Jkp ^ 2+B2cos{t '- r) p'dp'd0\ 



(2-69) 
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which conforms to the definition of a Bessel function. The same integration is 



performed in the second equation of (2-60), and the same Bessel function results. 
Finally, 



N„ = J 0 COS 6 COS 0 /r3 2 



2 



J v kayR 2 + B 2 
kayR 2 + Y 2 



J 



N c = -J 0 sin o Jia 2 



M ka\R 2 + B 2 

2 — 

kayR 2 + B 2 



(2-70) 



where J ] is a Bessel function of the first order. The far-field equations (2-57) 
become 



— cos 6 cos o na 2 



r 



i-jkr 



2jJ kaiR 2 + B 2 



kayR 2 + B 2 



r , Jkri e 2 

£ a = 1 sino Jia 2 

4/r r 



2 jJkayR 2 +B 2 



ka^R 2 + B 2 



(2-71) 



where 



R = u s -u = sin0 s cosp s - sin 0 co so 
B = u s -u = sin d s sin o s - sin d sino, (2-72) 



and 
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\fl 2 + B 2 = ^sin 2 o + sin 2 d s -2sin# s sin #cos(q - o s ). (2-73) 



Using equations (2-71) through (2-73) are computationally less 
intensive than a numerical integration. These equations were employed for both 
the far field pattern and gain computations. 

5. Near Field of the Antenna 

In the previous section, the far-field of a circular aperture antenna was 
determined. When a radome is located in the near-field, the integrals cannot be 
reduced to closed form and they must be evaluated numerically. The near field is 
generally considered to be the region surrounding the antenna where r << A , or 
kr « 1 [Ref. 10:pp. 106]. 

The complete expressions for the near-field of a circular aperture 
antenna excited by a current polarized in the x-direction are 



E» 

Ey 

E z 



j] + (* - e Jkr dx'dy' 



zJa 

*Kk 

-Jn 

4rtk 






-y)G 2 J H y jkr dx'dy' 

- z'"jG 2 J H ] e~ Jkr dx'dy\ 



(2-74) 



where 



Ar 2 r 2 - 1 - Jkr 
r 3 

3 + 3 jkr - k 2 r 2 



(2-75) 
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and 



r = - »' ) 2 + (y - y'f + (z - z'f . ( 2 - 76 ) 

The primed coordinates designate source points and the unprimed an observation 

point on the radome as in Figure 2.10. [Ref. 11] S' is the area of the antenna 
aperture. Since the radome is rotationally symmetric, the components \E H E y E z ) 

can be described in terms of cylindrical coordinates using the transformations 



\ 




Figure 2.10 Rectangular Coordinates for Circular Aperture Antenna 

= p'COS0' 

y = p'sinp' 

dx'dy' = p'dp'dy . (2-77) 

For an antenna scanned to (# 4 ,o s j the current distribution is 
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= J 0 [p') e 



jkx'ZQSQ'-® sin + y'sini o' - o sine 



( 2 - 78 ) 



Applying (2-75), 2-76) and (2-77) to equations (2-74) results in 



2 k a 
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0-0 p' = 0 



6, +(pCOS0-p'COS0') 2 G 2 
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2k a 
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y *rtk 
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[(pcosp-p'coso'Jjpsinp-p'sino'i^ p'dp'do' 



2k a 



E z =^ J J (p')e J * /, ' C0S ( 4 '* 1 )sine ' e~ Jkr x 



o =0 p =0 



( 2 - 79 ) 



|jp cos o - p' cos 0 '\jZG 2 p'dp'do'. 



where 



r = 




( 2 - 80 ) 



Since the expression for the excitation vector requires spherical components, 
[E H E y E z j are converted using a coordinate transformation 



34 



E r = E h sin dcos o + E y sin #sin o + E z cos 9 
E h = E h cos hcoso + E y cos 8 sin o - E z sin 0 
E 0 = -E h sin o + E y cos o. 



( 2 - 81 ) 



These are the equations programmed to compute the E-field impinging 
on the radome. The magnitude of the current distribution is allowed to vary as a 
function of p' . Various antenna sidelobe levels can be modeled by simply 
changing the function J 0 [p'\ 
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III. PROGRAMS 



The mathematical formulas developed in the previous chapters were 
computer programmed using FORTRAN. The pattern and gain calculations are 
done in two separate programs: RADOME.F and GAIN.F, both of which run on a 
Sun SPARCstation™ under UNIX. A brief description of each follows. 



A. DESCRIPTION OF THE PROGRAM RADOME.F 



This program, developed by Francis [Ref. 2], determines the series expansion 
coefficients of the current J s induced on the surface of the radome. It also 
calculates the total electric field surrounding the system, and generates radiation 
patterns using MATLAB™. Originally called LDBOR.F, (an abbreviation of 

"Loaded Bodies of Revolution") in Ref. 2, it has been modified to take advantage 
of symmetry, and also allows the antenna mainbeam to scan in both 0 and o . 



Figure 3.1 is a flow chart of the program. To run the program, data is entered 
interactively, except for a data file "gaus/gaus(ITH) M containing the integration 
constants for Gauss Quadrature, which is read in the program. All integrals are 
evaluated using this technique (see Appendix A.). Table I is a sample input. 

After all data is entered, the coordinate points of the surface are computed 
and stored in the arrays RH(= p) and ZH(= Z). The surface impedance of each 
segment is stored in ZLO. The elements of the MM impedance matrix are 
computed in subroutine ZMAT, and the excitation vector elements in GENEX. 



These provide 



■MM 



and 



U 



in equation (2-17). Subroutines DECOMP and 



SOLVE solve the matrix equation using Gaussian elimination. Subroutine 
PLANE calculates the measurement matrices, which are subsequently used to 
compute the scattered field from the radome. 
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Two antenna radiation subroutines are included. CIRCRTP calculates the 
electric field components and £, according to equations (2-81). It is used 

by subroutine GENEX to obtain the excitation vector. The second antenna 
subroutine ANTFF assumes that the observation point is in the far- field. This is 
used in the main program to compute the antenna contribution to the radiation 
pattern E a in (2-3). 

Output includes an ASCII file "outldbor" which lists all the input parameters, 
radome geometry and impedance, and the pattern magnitude. The current 
coefficients are written on a file "curcoefsdat” so that further patterns can be 
generated without recomputing the coefficients. Finally, MATLAB™ formatted 
files are generated for plotting the Q and o components of the electric field. 
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TABLE I. SAMPLE INPUT RADOME.F 



ENTER A LETTER TO INDIC ATE WHICH ITERATION THIS IS A 
SELECT BOR GEOMETRY BY NUMBER 
NUMBER BOR 

1 OGIVE 1 

2 SPHERE 

3 ('ONE 

4 DISK 

5 PARABOLA 

ENTER SURFACE CURVATURE (wavelengths) 7.5 

ENTER ZPRIME, WHERE CURVATURE STARTS (wavelengths) 0 
ENTER BASE RADIUS (wavelengths) 3 

ENTER FILENAMES gaus### 

ENTER THE FILENAME IN T (NT) gaus2 

ENTER THE FILENAME IN PHI (NPHI) gaus48 

ENTER THE FILENAME IN X AND Y gaus20 

ENTER THE PLOTTING INC REMENT IN DEGREES 1 

ENTER THE HIGHEST MODE 20 

ENTER PHI (observation) IN DEGREES 90 

ENTER SCAN ANGLES IN DEGREES: THETA, PHI 30,90 

ENTER COMPLEX IMPEDANCE. OHMS: (Real, Imag) (0,-1700) 

ENTER THE ANTENNA RADIUS (wavelengths) 1 
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Figure 3.1 RADOAfE.F Flow Diagram 
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Figure 3.1 RADOXIE.F Flow Diagram ( Continued) 
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Figure 3.1 RADOME.F Flow Diagram (Continued) 
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Figure 3.1 RADOME.F Flow Diagram (Continued) 
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Figure 3.1 RADOME.F Flow Diagram (Continued) 



44 




B. GAIN.F 



This program determines the gain of the system by solving (2-1 ). In order to 
perform the integration, the current coefficients from RADOME.F must be 
provided so that the radome scattered field can be computed. Input consists of 
specifying an iteration letter, and selecting the number of divisions in the 6 and 0 
integration. Input data is then read from a file "curcoefsdat" that is appended 
with the iteration letter chosen. 

Figure 3.2 is a flow chart of the program, and Table II is a sample input. 
GAIN.F contains several of the same subroutines that occur in RADOME.F. They 
include PLANE. CIRCRTP, and ANTFF. 

TABLE II. SAMPLE INPUT GAIN.F 



ENTER A LETTER TO INDICATE WHICH ITERATION THIS IS A 

ENTER NUMBER OF DIVISIONS IN THETA DIRECTION, 

NDIVT = ? 2 

ENTER NUMBER OF DIVISIONS IN PHI DIRECTION, NDIVP =? 3 



( Begin ) 

nr 



Assign Parameters 



i 



Read ITH,S1 LECTION, BASE 
NP,MT,MP.N.NR( )WS, 
MODES,NBLOCK,THS,PHS, 
ARAD, PHI. IMP.ZP,RS,RH, 
ZH,ZLO from "curcoefsdatdTH)" 



data was written to M curcoefsdat(ITH) M 
from RADOME.F 



! RoadCT(L) 



L 



CT(L) are elements of current coefficient 
vector [I] written to curcoefsdat 



Read XT(I), AT((). X(I), A(I), XP(I), AP(I) 



A 



7 



A(I) are weights, and X(I) are adscissas 
which are read from gaus/gaus# 



Enter NDIVT, NDf VP 



7 



NDIVT, and NDIVP are the number of 
divisions in 0 and 0 respectively. This 
allows for easy variation of the number 
of integration steps without recomputing 
weights and abscissas for Gauss 
integration. 



Calculate SO ), b integration angles 
and Q(I), 0 integration angles 



I 



Set SUM, SUMNONE. UMAX. UMAXNONE, 
EMAX, EM AXNC )NE to 0 




Figure 3.2 GAIN.F Flow Diagram 
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o 



Figure 3.2 GAIN.F Flow Diagram ( Continued) 
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Figure 3.2 GA/N.F Flow Diagram ( Continued) 
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IV. RESULTS AND CONCLUSIONS 



A. VALIDATION AND CALCULATED RESULTS 

Part of the continued development of this research is validation of the 
programs. Test cases were run to see how the programs performed against 
known results. Test radomes ranged from nearly nonexistent disks located at 
large distances, to perfectly conducting parabolas located in the near-field. E- 
plane and H-plane patterns defined in Figure 4.1 were plotted to determine the 



E-plane H-plane 





Figure 4.1 Orientation of E-plane and H-plane 



effects of the radome. The antenna pattern without the radome was plotted 
(dashed line) along with the pattern of the complete system (solid line) to 
illustrate the effects of radome scattering. 
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1. Code Validation: Scattering from a Small Disk 

The pattern of an antenna radiating in the presence of a negligible -small 
circular-disk located at a large distance was computed. Since the radome 
scattered field is negligible, the predicted output is that of a uniform!) -excited 
circular-aperture antenna. Figure 4.2 shows the plots of the E-plane pattern for 
three different radius antennas (.5 X. 1.0 X and 2.0 X). The plots appear as single 
lines with no scattering. The gain as computed by GAIN.F. was compared to the 
theoretical gain of a circular aperture with a uniform current distribution 



where 3 is the aperture radius [Ref. 10: pp. 483 J. Theoretical gains for the .5 X. 
1.0 X and 2.0 X radius apertures are 9.94. 15.96. and 21.98 dB respectively. In 
each case the numerical gain is within 1.5 dB of the theoretical gain. 

2. Code Validation: Perfectly Conducting Paraboloid 

Programs have been written that use the MM solution for bodies of 
revolution developed by Mautz and Harrington to predict the patterns and gains 
of reflector antennas [Ref. 7|. Assuming that the radome is a perfect! v 
conducting parabola, the output of RADOME. F and GAIN.F can be compared 
with the output of these existing programs. 

Figure 4.3 is a test case where the radome is a parabola with surface 
impedance Z, = 0 + J0 Q . The dimensions of the radome matched those of a 
reflector which was analyzed by a separate computer program that calculates 
both pattern, and gain. The pattern and gain for the perfectly conducting radome 



2 




(4-1) 
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60 1 — f [ ' — I 

0 50 100 150 200 250 300 350 



Theta Angle (deg) 

Figure 4.2 Scattering from a Negligibly Small Disk 
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computed by RADOME.F and GAIN.F are identical to that computed for the 
reflector. 





3. Dielectric Radomes 

To determine the effect of various radome materials on antenna 
performance, radomes with several different surface impedances were analyzed. 
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Surface impedances were calculated using (2-32) for a range of e r and tan 8 
typical of radome materials. The real and imaginary components of Z L as a 
function of e r and tan 8 are plotted in Figure 4.4. As shown in the figure, the 
reactance is essentially independent of tan<5 over the range 0<tan<5<0.1, 
and the real component of Z L is strongly dependent on tan <5, particularly in the 
region of low e r . 

A test radome consisting of an ogive with dimensions given in Table I 
and parameters (l x = l/28 and tan 5 = 0.01 was analyzed.* The patterns were 

calculated for dielectric constants ranging 2 < e r < 10 , and the results are 
summarized in Table III. The E and H-plane patterns of an ogive radome with 
Z L = 34- J1700 Q are shown in Figures 4.5 and 4.6. The impedance 
corresponds to a radome made out of a material of e r =2.0 and tan <5 = 0.01. 
The scattering is symmetric with respect to 0 for the unscanned antenna, and the 

backlobe level is -45 dB below the peak. The scanned case has higher sidelobes. 
as expected, with a lobe of -30 dB at an angle 270°. The gain loss due to the 
radome for the unscanned case is approximately 0.1 dB, and for the scanned case 

0.2 dB. The discontinuity in the H-plane patterns at 0 = 90° is due to the 
antenna model. When <j> = 90° the field is § polarized, and according to (2-71). 

E a has no cos# factor to force the pattern to 0 at # = 90°. Note that the 

radome ohmic loss has not been included. 

Increasing the dielectric constant results in a constant-shape scattering 
pattern with increasing magnitude. Figure 4.7 is the E-field pattern for e r =10, 
which corresponds to some ceramics. It is apparent that the scattering by the 



* These dimensions are not intended to match any particular existing 
radome. 
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radome in the rear hemisphere is significant, and the gain loss is nearly 4 dB. 
Table 111 shows the gains with and without the radome for 2 < e r < 1 0. Figure 
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Figure 4.4 Surface Impedance Z L as a Function of e r and tan 8 
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4.8 is a plot of the gain loss relative to no radome as a runction of e r . For 
reference the Fresnel reflection coefficient for a planar interface is also shown. 
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TABLE III. GAIN LOSS FOR AN OGIVE RADOME 



^ r 


( n ) 

(tan<5 = 0.01 J 


Gain (dB) 
WithAVithout 
Radome 

e s =0° (o, =0°) 


Gain (dB) 
WithAVithout 
Radome 

9 S = 30° (c) s = 0°) 


2.0 


0-j 1700 


22.21/22.32 


21.61/21.78 


2.0 


34-j 1700 


22.21/22.32 


21.61/21.78 


4.0 


7.5-/560 


21.42/22.32 


20.90/21.78 


5.0 


5.3-j420 


20.91/22.32 


20.46/21.78 


7.0 


3.3-j280 


19.95/22.32 


19.47/21.78 


8.0 


2.7-j240 


19.51/22.32 


18.96/21.78 


9.0 


2.4-j210 


19.04/22.32 


18.47/21.78 


10.0 


2.1-jl90 


18.59/22.32 


18.00/21.78 
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Figure 4.5 E- Plane Pattern of Antenna with an Ogive Radome 

Z ( = 34 - j\ 700 Q 
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Theta Angle (deg) 

Figure 4.6 H-Plane Pattern of Antenna with an Ogive Radome 

Z L = 34 - j\ 700 Q 



58 





Figure 4.7 E-Plane Pattern of Antenna with an Ogive Radome 

Z t = 2. 1 - j\ 87 Q 
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The gain loss appears to be linear with e r . This is probably because the base of 
the ogive is open, which allows reflected waves to exit the radome. If a back 
surface was present, these waves would be reflected again causing additional 
gain loss. 

a. Conical Radome 

To determine the effect that radome shape has on scattering, a 
conical radome was analyzed. The length, width, and impedance of the cone was 
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Figure 4.8 Gainloss Relative for an Ogive Radome as a Function of e r IE- 

Plane Scan) 



the same as the ogive radome discussed previously. A conical radome presents a 
relatively flat surface that could possibly cause more focused scattering than a 
curved shape like an ogive. Ray tracing of Figure 4.9 shows that the reflections 
from the radome will add coherently for the cone, and non-coherently for the 
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ogive. However, this assumes that the radoine is in the far field of the antenna. 
The E-plane pattern is shown in Figure 4.10. 



x 




x 




Figure 4.9 Scattering from Conical Radome vs. Ogive 

A comparison of Figures 4.5 and 4.10 demonstrates that there is 
even less scattering from the conical radome in the rear hemisphere. For the 
scanned case, the conical radome has a distinct lobe at 9 = 280° due to the 
specular reflection of the main lobe from one side of the radome. This is referred 
to as a reflection lobe. 
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Figure 4.10 E-Pane Pattern of Antenna with a Conical 
RadomeZ t = 34 - jl 700 Q 
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B. CONCLUSIONS 

This research has resulted in two computer codes that can be used to predict 
the pattern degradation and gain loss due to an axially symmetric radome in the 
near field of a linearly polarized circular aperture. The computed patterns show 
that the gain loss is slightly greater than the reflection loss for a planar interface 
with the same dielectric constant as the radome. This is approximately true for 
low-loss radomes. For lossy radomes the impedance has a very strong affect on 
scattering, while radome shape has less effect. Scanning the antenna results in an 
asymmetric scattering pattern and generally higher sidelobe levels. 

The computed gain compares favorably with values computed by closed form 
solutions derived from theory, and the gain of the system has the expected 
decrease as the level of scattering increases. For the specific case of the ogive 
radome, there is a linear dependence of the gain loss to the relative dielectric of 
the radome. Note that the gain calculations presented here do not include ohmic 
loss inside of the radome. 

Improvements in the original code include arbitrary aperture illumination and 
incorporating mode symmetry. Computational run time is reduced by nearly half 
when symmetry is used, and significantly increases the applicability of this 
method of analysis. 

In order to apply this analysis to multilayer radomes, an equivalent 
surface impedance must be determined. The program can be modified to account 
for radomes constructed of nonuniform materials by assigning different 
impedances for each segment of the plane curve generating the BOR. 
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APPENDIX A. GAUSS QUADRATURE 



Computing the MM elements of 



U 



and 



(2-19 and 2-28) as well as the 



gain requires the numerical evaluation of integrals. MM requires numerous 
integration in order to solve for the current coefficients l' nj and l° nJ of (2-31). Not 

only is numerical analysis necessary to perform the integrals in the MM. but also 
the matrix analysis of (2-3 1 ). The large number of separate integrals necessitates 
an integration technique which provides accuracy in a few steps to avoid 
excessive run time. Gauss Quadrature is a method of numerical integration which 
employs unequally spaced intervals . and provides accuracy with a few points. 
[Ref. 4:pp. 154-157] 

Gauss Quadrature approximation of the form 



b 

I = \f(x)dx 



a 



* ^ 



(A-l) 



by the sum of m terms 



m 

I = *7* I U>Jl « k ). 

k = 1 



(A-2) 



The coefficients UJ k are weights, and the H k s are abscissas determined from 

b + a b-a K . » 

+ — ( A - 3 > 



where % k ' s are the m number of zeros of the m'th degree Legendre Polynomials. 



64 



APPENDIX B. 



RADOME.F LISTING 



The material on the following pages is a listing of the program RADOME.F 
described in chapter III. 
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c************************************************************************ 
C ROGRAM radome.f V.3 

C DATE : 23 January 1992 

C LATEST REVISION: 31 December 1992 

C PROGRAMMERS D. JENN, R. FRANCIS, K. KLOPP 

C 

C NOTES: 

C 1. ARBITRARY CIRCULAR APERTURE ILLUMINATION SPECIFIED IN 
C FUNCTION TAPER. 

C 2. NEAR FIELD COMPONENTS INCLUDED. (COMPUTED IN SUBROUTINE GENEX . ) 

C 3. ANTENNA IS AT THE BOR COORDINATE SYSTEM ORIGIN. 

C 4. LINEARLY (X-POLARIZED) ANTENNA. 

C 5. CLOSED FORM FAR FIELD ANTENNA OPTION ADDED. 

C 6. FLAGS USED THROUGHOUT: (IMP IS SET AUTOMATICALLY) 

C IMP=0 perfect conductor radome (for test purposes) 

C IPRINT=0 print pattern points to unit 8 

C ISEG=0 print the generating curve points to unit 8 

C ICALC=0 compute current coefficients & pattern 

C #0 read current coefficients from disc file given 

C by ’ CURCO and compute the pattern. 

C IFF=0 closed form far field antenna pattern is used to 

C compute ETF and EPF (SUBROUTINE ANTFF ) 

C #0 far field computed using CIRCRTP at distance ROBS 

C ISYM=0 use mode symmetry 

C #0 turn off mode symmetry (for testing) 

C 7. THE FACTOR ETA (=377 OHMS) IS OMITTED THROUGHOUT. 

C 

C BEGIN MAIN PROGRAM********************************************* 
CHARACTER ITH 

CHARACTER*9 OUT , ETSCATTER , EPSC ATTER , CURCO 
CHARACTER*12 CC 
CHARACTER*6 AN , ETFEED , EPFEED 
CHARACTER*8 GNT , GNPHI , CGP , CPOL , XPOL 
CHARACTER*14 TPTS , PPTS , PHIPTS 

COMPLEX EP, ET,Z( 100000 ) ,R( 1600) ,B(800) ,C (800) ,U 
COMPLEX UC , ET 1 , EP 1 , ET2 , EP2 , EC , EX , ZLO ( 400 ) , ZL ( 2400 ) 

COMPLEX CIRCR, CIRCT , CIRCP , ETF , EPF 

COMPLEX EXP1 , EXP2 ,CONJG , CEXP , CMPLX , CT(30000) , IMP , JK 
DIMENSION RH(400) , ZH (400 ) , XT(4) , AT(4) , IPS(800) ,XP( 100) , AP( 100) 
DIMENSION A(lOO) , X( 100) , EXP (500) , ANG (500) , ECP (500) 

DIMENSION ECV (500) , EXV(500) ,PHC(500) , PHX(500) 

DIMENSION ETSCAT(500) , EPSC AT (500) , ETHF(500) , EPHF(500) 

INTEGER CNPHI, SELECTION 

DATA IPRINT/O/, ICALC/l/, IFF/O/ , ISEG/O/ , ISYM/O/ 

DATA START, STOP/O. , 360 ./, PI/3 . 1415926/ 

Rad=PI/ 180 . 

ECX=0 . 

BK=2 . *PI 
U=(0. ,1.) 

U0=(0. ,0. ) 

UC=-U/4 . /PI 
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JK=(0. 0,6. 283185307) 

enter a letter to indicate which iteration this is ****«»*«.*** 

C***»* (ALL DATA FILES ARE APPENDED WITH THIS LETTER) *********** 

C***«* (.m FILES ARE FOR GENERATING PLOTS IN MATLAB) *********** 

WRITE (6,*) ’ENTER A LETTER TO INDICATE WHICH ITERATION THIS IS’ 
READ(5, *)ITH 
OUT= ’ outldbor ’ //ITH 
AN = ’ ang ’//ITH/ / ’ .m’ 

ETFEED=’etf ’//ITH//’ ,m’ 

EPFEED= ’ epf ’ //ITH// ’ . m ’ 

ETSCATTER= ’ etscat ’//ITH// ’ .m’ 

EPSCATTER= ’ epscat ’//ITH/ / ’ .m’ 

CURCO=’asciidat '//ITH 



CC= ’ curcoef sdat ’ //ITH 
CP0L= ’ etpol ’//ITH// ’ .m’ 

XP0L= ’ eppol ’ //ITH/ / ’ .m’ 

IF(ICALC.Eq.O) THEN 

WRITE (6,*) ’ENTER THE ANTENNA RADIUS (wavelengths)’ 

READ (5 , *) ARAD 

(^♦♦♦♦♦SELECT THE GEOMETRY OF THE BOR************************************* 



WRITE ( 6 , * ) 
WRITE ( 6 , * ) 
WRITE(6,*) 
WRITE(6,*) 
WRITE (6 , *) 
WRITE(6,*) 
WRITE(6 ,*) 



SELECT BOR GEOMETRY BY NUMBER. 



NUMBER BOR’ 

1 OGIVE’ 

2 SPHERE’ 

3 CONE’ 

4 DISK’ 

5 PARABOLA’ 

READ (5,*) SELECTION 

IF (SELECTION. EQ. 1) THEN 

WRITE(6,*) ’CALLING SUBROUTINE FOR OGIVE’ 

CALL OGIVE ( NP ,ZH , RH , BASE , RS ,ZP) 

ELSEIF ( SELECTION . EQ . 2)THEN 

WRITE(6,«) ’CALLING SUBROUTINE FOR SPHERE’ 

CALL TESTSPHERE(NP , ZH , RH , BASE , RS , ZP) 

ELSEIF ( SELECTION . EQ . 3 )THEN 

WRITE(6,*) ’CALLING SUBROUTINE FOR CONE’ 

CALL CONE(NP , ZH , RH , BASE , RS ,ZP) 

ELSEIF( SELECTION . EQ . 4 )THEN 

WRITE(6,*) ’CALLING SUBROUTINE FOR DISK’ 

CALL DISK(NP , ZH ,RH , BASE, RS , ZP) 

ELSEIF ( SELECTION .EQ. 5) THEN 

WRITE(6,*) ’CALLING SUBROUTINE FOR PARABOLA’ 
CALL PARAB ( NP , ZH , RH , DM , FOD , ZP ) 

BASE=DM/2 . 

ENDIF 

IF (NP . GT. 400)THEN 

WRITE( 6,*) ’MAXIMUM NUMBER OF POINTS(NP) IS 400’ 
GOTO 999 
ENDIF 
ENDIF 



C* *********************************************** ************ 
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WRITE(6,*) ’ENTER THE FILENAMES gaus###’ 

WRITE (6 , *) ’ ’ 

WRITE(6,*) ’ENTER THE FILENAME IN T (NT)’ 

READ(5,*)GNT 

WRITE(6,*) ’ENTER THE FILENAME IN PHI (NPHI) ’ 

READ (5 , *)GNPHI 

WRITE(6,*) ’ENTER THE FILENAME IN X AND Y’ 

READ(S , *)CGP 

C OPEN THE FILES FOR THE gaus/gaus### 

TPTS= ’ gaus/ ’ / / GNT 

PPTS= ’ gaus/ ’ / / GNPHI 

PHIPTS=’gaus/’//CGP 

OPEN ( 1 , FILE=TPTS , STATUS= ’ OLD ’ ) 

OPEN ( 2 , FILE=PPTS , STATUS^ ’ OLD ’ ) 

OPEN (4 , FILE=PHIPTS , STATUS= ’OLD ’ ) 

READ(1 ,*)NT 

IF(NT.GT.4)THEN 

WRITE(6,*) ’MAXIMUM NUMBER OF POINTS(NT) IS 4’ 

GOTO 999 
END IF 

READ(2 , *) NPHI 

IF(NPHI.GT. 200) THEN 

WRITE(6,*) 'MAXIMUM NUMBER OF POINTS(NPHI) IS 200’ 
GOTO 999 
ENDIF 

READ(4 , *)CNPHI 

IF(CNPHI . GT. 100)THEN 

WRITE(6 ,*) ’MAXIMUM NUMBER OF POINTS(CNPHI ) IS 100’ 
GOTO 999 
ENDIF 

C LOAD THE WEIGHTS AND ABSCISSAS IN THE VECTORS. 

DO 1 M= 1 , NT 

READ(1,*,END=1)XT(M) ,AT(M) 

1 CONTINUE 

DO 2 M=1 , NPHI 

READ ( 2 , * , END =2 ) X (M ) ,A(M) 

2 CONTINUE 

DO 4 M=1 , CNPHI 

READ(4,*,END=4)XP(M) ,AP(M) 

4 CONTINUE 
CLOSE(l) 

CL0SE(2) 

CL0SE(4) 

WRITE (6,*) ’ENTER THE PLOTTING INCREMENT IN DEGREES’ 

READ (5 , *)DT 
IF(ICALC.EQ.O) THEN 
WRITE(6,*) ’ENTER HIGHEST MODE’ 

READ (5,*) MODES 

NBL0CK=2*M0DES+1 

NR0W=NBL0CK*(2*NP-3) 

C CHECK FOR BOUNDS VIOLATION 
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IF(NROW.GT. 30000) THEN 

WRITE(6 , * ) ’BOUNDS VIOLATION FOR ARRAY CT(.)’ 

GOTO 999 
ENDIF 
ENDIF 

WRITE(6,*) ’ENTER PHI (observation) IN DEGREES’ 

READ(S , *)P 
PHI=P*RAD 

IF(ICALC.EQ.O) THEN 

WRITE(6,*) ’ENTER THE SCAN ANGLES IN DEGREES: THETA , PHI ’ 

READ(5 , *) THD.PHD 

THS=THD*RAD 

PHS=PHD*RAD 

WRITE(6,*) ’ENTER COMPLEX IMPEDANCE, OHMS: (Real , Imag) ’ 

READ (5 , *) IMP 
ELSE 

C**********HEAD CURRENT COEFFICIENTS IF ICALC80************************ 
WRITE(6 , *) ’ENTER LETTER CODE FOR FILE curcoefsdat’ 

READ(S , *) ITH 
CC= ’ curcoefsdat ’ //ITH 
0PEN(31 ,FILE=CC) 

WRITE(6 ,*) ’FILE OPENED IS \CC 

READ (31,*) ITH, SELECTION, BASE, NP,MT,MP,N 

READ (31,*) NROW , MODES , NBLOCK , THS , PHS , ARAD , PHIO 

READ (31 , *) IMP ,ZP ,RS 

READ (31 , *) (RH(L) ,L=1,NP) 

READ (31 , *) (ZH(L),L=1,NP) 

READ (31 ,*) (ZLO(L) ,L=1 ,NP) 

READ(31 , *) (CT(L),L=1, NROW ) 

CL0SE(31) 

URITE(6 , *) ’RUN CODE READ FROM DISC : ’,ITH 
ENDIF 

c**********ready to calculate the pattern***************************** 

IT=INT( (ST0P-START)/DT)+1 

MHI=M0DES+1 

OPEN ( 8 , FILE=OUT ) 

WRITE( 8 , 8000 ) 2 . *BASE , NP , MODES , RS , ZP 
8000 FORMAT(//, ’*** BOR RADIATION PATTERN FOR A CIRCULAR’ 

*’ DISC RADIATOR USING GENEX ***’, 

*// , 2X , ’ BOR DIAMETER (WAVELENGTHS) 3 ’ ,F5.2,/,2X, 

* ’NUMBER OF GENERATING POINTS (NP)= * , 14,/, 2X , 

* ’HIGHEST MODE NUMBER USED (MODES) 3 ’, 13 , 

*//,2 X, ’SURFACE RADIUS’ ,FS. 2,/, 2X, ’ZPRIME’ ,F5.2) 

WRITE (8, 30) NT, NPHI , CNPHI 
30 F0RMAT(/,12X, ’ NT NPHI CNPHI’, 

* / , 10X , 15 , 2X , 15 , 2X , 15) 

IF(ISEG.EQ.O) WRITE(8 , 1300) 

1300 F0RMAT(/,10X, ’INDEX’ ,8X, ’Z(I) ’ , 10X , ’ RHO(I) ' ,12X, ’SURF IMPED’ ) 
MP=NP-1 
MT=MP-1 
N=MT+MP 
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DO 52 1=1, NP 

IF(ABS(ZH(I)).LT. .001) ZH(I)=0. 

IF(ABS(RH(I)) .LT. .001) RH(I)=0. 

ZHB=ZH(I)/BK 

RHB=RH(I)/BK 

C ASSIGN SURFACE IMPEDANCE AT THIS POINT. THE SURFACE IMPEDANCE OF SEGMENT 
C I IS ZLO(I) 

IF(ICALC.EQ.O) Z10(I)=IMP/ ( 120 . *PI ) 

IF(ISEG.Eq.O) WRITE(8 , 8004) I , ZHB ,RHB , IMP 
52 CONTINUE 

8004 FORMAT( 1 IX , 14 , 4X , F8 . 3 , 8X , F8 . 3 , 6X , 2F8 . 2) 

Qm********** ********************************* **************************** 

C BIG LOOP *>***************************** 

c 

C MODE LOOP TO CALCULATE THE CURRENT COEFFICIENTS. POS AND NEG MODES * 
:ONE IN THE SAME ITERATION OF THE LOOP 



C * v 

C 400 CONTINUE < **************************** 



C*********************** * ****** ************* ***£ LOAD , ZMAT , GENEX , DECOMP , SOLVE 
IF(ICALC.EQ.O) THEN 

IF(CABS(IMP) .NE.O) CALL ZLOAD(NP , RH , ZH , ZLO , ZL) 

DO 400 M=1 , MHI 
NM=M-1 

CALL ZMAT(NM , NM , NP , NPHI , NT , RH , ZH , X , A , XT , AT , Z ) 



C**** *********************************************** ****POSITIVE MODE +n 
C I t,t t,phi 

C I Z Z 

C | +n '+n 

C ZMAT = 

C 
C 
C 

C*********************************************************************** 

C MODE SYMMETRY IN n FOR ZMAT EXISTS BUT IS NOT USED IN THIS VERSION. 

C MODE SYMMETRY IN I USED DIRECTLY FOR PHIS OF 0,90,180, AND 270. 

O********************************************************************** 



phi,t phi, phi 
Z Z 

+n +n 



t ,t 

Z 

-n 



t ,phi 

Z 

-n 



phi,t phi, phi 
Z Z 



Z 

+n 



t ,phi 
-Z 
+n 



phi,t phi, phi 

-Z Z 

+n +n 

O*** ************************************************ CURRENT COEFFICIENTS 

C -1 

C 

c 

C I +n I ■ I I +n 

C 
C 

C CT = 
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0 



0 



c 
c 
c 
c 
c 
c 

C************************************************ * ** POSITIVE MODE +n 



-n 



i. 



.1. 




-n 



IF(CABS(IMP) .NE.O) CALL ZTOT(MT,MP,ZL,Z) 
CALL GENEX (NM , NP , NT , NPHI , CNPHI , XT , AT , X , A , 

* xp,ap,ths,phs,arad,rh,zh,b) 

CALL DECOMP(N , IPS , Z) 

CALL SOLVE(N,IPS,Z,B,C) 

NTOPl=MODES-NM 
NT0P2=NBL0CK-(NT0P1 +1 ) 

NS2=NT0P1*N 

NS1=NT0P2*N 

NROW=NBLOCK*N 



C*************************************************************** 

C STORE CURRENT COEFFICIENTS IN ONE LONG VECTOR CT(.). MOST NEGATIVE 
C MODE IS AT THE TOP; MOST POSITIVE AT THE BOTTOM. UNIT 30 IS FORMATED 
C OUTPUT; UNIT 31 IS FREE FORMATED (USED BY GAIN PROGRAM). 

C************************* ************************************** 

OPEN ( 30 , FILE=CURCO ) 

WRITE (30, 5027) ’Mode( ’ ,NM, ’ ) = ’ , NM , ’ NM = ’ ,NM 
DO 401 L=1 , N 

WRITE(30 , 5026) ’CT(’,L+ NS1,’) = ’ ,C(L) 

401 CT(L+NS1)=C(L) 

5027 FORMAT (/A, I,A,I/A,I/) 

5026 FORMAT (5X ,A,I4,A,F10.4, ’+ ( ’ ,F10 . 4, ’ )*i ’ ) 

C***************END OF POSITIVE MODE*******BEGIN NEGATIVE MODE -n 
C MODE SYMMETRY IS USED IF ISYM=0 
IF(NM.NE.O) THEN 
NMN=-NM 

C USE BRUTE FORCE IF ISYM # 0 
IF(ISYM.NE.O) THEN 

CALL GENEX(NMN , NP , NT , NPHI , CNPHI , XT , AT ,X , A , 

* XP.AP.THS.PHS.ARAD.RH.ZH.B) 

CALL ZMAT(NMN ,NMN ,NP ,NPHI ,NT, RH , ZH , X , A, XT, AT, Z) 

IF(CABS(IMP) .NE.O) CALL ZTOT(MT,MP,ZL,Z) 

CALL DECOMP(N.IPS.Z) 

CALL SOLVE(N,IPS,Z,B,C) 

ENDIF 

WRITE(30,5027) ’Mode( ' ,NMN, ’ ) = '.NMN.’NM = * ,NM 

DO 402 L=1 ,MT 

CT(L+NS2)=+C(L) 

WRITE(30,5026) ’CT( ’ .L+NS2, 1 ) = ’,CT(L+NS2) 

402 CONTINUE 

DO 403 L=MT+1,N 
IF(ISYM.NE.O) CT(L+NS2)-+C(L) 

IF(ISYM.EQ.O) CT(L+NS2)=-C(L) 

WRITE(30 , 5026) ’CT( ’ .L+NS2 , ’ ) = ’ ,CT(L+NS2) 
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403 CONTINUE 

C******************************************END OF NEGATIVE MODE -n 
ENDIF 

400 CONTINUE 
CL0SE(30) 

q****** ***************************************** ************************ 

C WRITE THE VECTOR OF CURRENT COEFFICIENTS ON DISC FOR GAIN CALC. 

C (FREE FORMAT TO UNIT 31.) 

£***************************♦******************************************* 
OPEN (31 ,FILE=CC) 

WRITE(3 1 , * ) ITH , SELECTION , BASE , NP , MT , MP , N 
WRITE (31,*) NROW , MODES , NBLOCK , THS , PHS , ARAD , PHI 
WRITE(31 , *) IMP , ZP ,RS , IFF 
WRITE(31 , * ) (RH(L) , L=1 ,NP) 

WRITE(31 , *) (ZH(L) ,L=1 ,NP) 

WRITE(31 , *) (ZLO(L) ,L=1 ,NP) 

WRITE(31 , *) (CT(L) ,L=1 ,NROW) 

WRITE(31 , *) GNT.GNPHI ,CGP 

WRITE(6,*) ’CURRENT COEFFICIENTS WRITTEN ON DISC’ 

CL0SE(31) 

ENDIF 

Q* ************************************** ********************************* 

C END OF BIG LOOP *>************************* 
C * v 

C*******< ************ *************** 

Q*** *************************************************************** ****** 

c 

C BEGIN PATTERN LOOP FROM START TO STOP IN INCREMENTS OF DT (ALL IN DEG) 

C 

DO 500 1=1, IT 

THETA=START + DT*FLOAT( 1-1 ) 

THX=THETA*RAD 

PHR=PHI 

IF(THETA.GT. 180. ) THEN 
PHR=PHI+PI 

THX= (360. -THETA )*RAD 
ENDIF 

ET1=(0 . ,0. ) 

EP1= (0 . ,0. ) 

ET2= (0 . ,0. ) 

EP2=(0 . ,0. ) 

DO 300 M=1 ,MHI 
NM=M- 1 

EXP1=CEXP(CMPLX(0. , FLOAT (NM) *PHR) ) 

EXP2=C0N JG (EXP 1 ) 

C**** ****************************************************** * PLANE 
CALL PLANE( NM , NH , NP , NT , RH , ZH , XT , AT , THX , R) 

NT0P1=M0DES-NM 
NT0P2=NBL0CK-(NT0P1+1 ) 

NS2=NT0P1*N 

NS1=NT0P2*N 
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DO 250 L=1,MT 

ET1=ET1+R(L)*CT(L+NS1)*EXP1 
EP1=EP1+R(L+N) *CT(L+NS1 ) *EXP1 
IF(NM.EQ.O) GOTO 250 
ET1=ET1+R(L)*CT(L+NS2)*EXP2 
EP1=EP1-R(L+N) *CT(L+NS2) *EXP2 
250 CONTINUE 

DO 260 L=1,MP 

ET2=ET2+R(L+MT ) *CT(L+NS1+MT) *EXP1 
EP2=EP2+R(L+MT+N)*CT(L+NS1+MT)*EXP1 
IF(NH.EQ.O) GO TO 260 
ET2=ET2-R(L+MT)*CT(L+NS2+MT)*EXP2 
EP2=EP2+R(L+MT+N)*CT(L+NS2+MT)*EXP2 
260 CONTINUE 
300 CONTINUE 

C FEED CONTRIBUTION IN THE FAR FIELD IS ETF , EPF 

C***» ********************************************* ***gESSJ 1 . 

C USE CLOSED FORM FEED EXPRESSION FROM SUBROUTINE ANTFF IF 
C IFF=0 ; OTHERWISE USE BRUTE FORCE EVALUATION FROM CIRCRTP 
C 

R0BS=1000 . 

IF(IFF.EQ.O) THEN 

CALL ANTFF (THX , PHR , THS , PHS , ARAD , ETF , EPF) 

ELSE 

RHB=ROBS*SIN (THX) 

ZHB=R0BS*C0S(THX) 

CALL CIRCRTP(CNPHI , XP ,AP, ARAD, THS, PHS, 

* _ PHR , RHB , ZHB , CIRCR , CIRCT , CIRCP ) 

C REMOVE THE 1/R DEPENDENCE BECAUSE EXP(-jkR)/R IS OMITTED IN 
C THE SCATTERED FIELDS ET AND EP 
ETF=CIRCT*ROBS 
EPF=CIRCP*ROBS 
ENDIF 

ETHF ( I ) =C ABS (ETF ) 

EPHF(I)=CABS(EPF) 

ET=(ET1+ET2)*UC 

EP=(EP1+EP2)*UC 

ETSCAT(I)=CABS(ET) 

EPSCAT( I) =CABS(EP) 

C TOTAL E-THETA AND E-PHI COMPONENTS 
EC=ET+ETF 
EX=EP+EPF 
ECV(I)=CABS(EC) 

EXV(I )=CABS(EX) 

ECR=REAL(EC) 

ECI=AIMAG(EC) 

EXR=REAL(EX) 

EXI=AIMAG(EX) 

PHC(I)=ATAN2(ECI , ECR+1 . e-10) /RAD 
PHX(I) =ATAN2(EXI ,EXR+1 . e- 10) /RAD 
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ANG(I)=THETA 

ECX=AMAX1 (ECX , ECV( I ) ,EXV(I)) 

500 CONTINUE 

WRITE(6 , *) ’MAX E VALUE=’,ECX 
WRITE( 8 , 103 ) P , THS/RAD , PHS/RAD , ARAD , ECX 
103 FORMAT(/, 10X, ’PHI OF RECEIVER (DEG)= ’ , F8 . 2, /, 10X , 

* ’ANTENNA SCAN: THETA (DEG) = ’ , F8 . 2 , / , 10X , 

* ’ PHI (DEG)=’ ,F8.2,/,10X, 

* ’ANTENNA RADIUS: ARAD= ’ ,F8 . 3,/ , 10X , 

* 'MAXIMUM FIELD VALUE (V/M)= ’ ,E15 . 5) 

IF (IFF . EQ . 0) THEN 

WRITE(8 ,113) 

ELSE 

WRITE(8 , 213) ROBS 
ENDIF 

113 F0RMAT(/,10X, ’CLOSED FORM FAR FIELD PATTERN FROM ANTFF IS USED’) 
213 F0RMAT(/,10X, ’FAR FIELD COMPUTED USING CIRCRTP, ROBS=’ ,F10.2) 

C STORE DATA FOR NORMALIZED CO- AND CROS-POLARIZED PATTERNS 
DO 600 1=1, IT 
ECP(I)=(ECV(I)/ECX)**2 
EXP(I)=(EXV (I)/ECX)**2 
ECP(I)=AMAX1 (ECP(I) , 1 .E-6) 

EXP(I)=AMAX1(EXP(I) .l.E-6) 

ECP(I)=10. *AL0G10(ECP(I ) ) 

EXP(I)=10 . *AL0G10(EXP(I ) ) 

600 CONTINUE 

IF(IPRINT.EQ.O) THEN 
WRITE(8 ,5015) 

5015 F0RMAT(//,7X, ’ANGLE’ , 15X, ’CO-POL’ ,25X, ’X-POL’ ,/,7X , 

1 ’ (DEG) ’ ,4X,2( ’(VOLTS) ’ ,4X, ’(DEG) ’ ,3X, ’ (DB-REL) ’ ,4X)) 

DO 9000 L=1 , IT 

WRITE(8,5016) ANG(L) ,ECV(L) ,PHC(L) ,ECP(L) ,EXV(L) ,PHX(L) 

1 ,EXP(L) 

5016 FORMAT ( 5X ,F6 . 2 , 3X ,2(F8 .4 , 3X , F7 . 2 ,3X ,F7 . 2 , 3X) ) 

9000 CONTINUE 

ENDIF 

0PEN(2,file=AN) 

0PEN(3 , f ile=CPOL) 

OPEN (4 ,f ile=XPOL) 

0PEN(7 , lile=ETSCATTER) 

OPEN (8 ,lile=EPSCATTER) 

OPEN (9 , f ile=ETFEED) 

OPEN ( 10 , f ile=EPFEED ) 

DO 9097 1=1, IT 
WRITE(2, 5019)1, ANG(I) 

WRITE(3, 5020)1, ECP(I) 

WRITE(4, 5021)1, EXP(I) 

WRITE(7, 5022)1, ETSCAT(I) 

WRITE(8, 5023)1, EPSCAT(I) 

WRITE(9, 5024)1, ETHF(I) 

9097 WRITE(10, 5025)1, EPHF(I) 
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5019 FORMAT ( ' ANG( ’ , I , ’ , F8 . 3 , ’ ; ’ ) 

5020 FORMAT ( ’ ECP ( ’ , I ,')=*, F8 . 3 , ’ ; ' ) 

5021 FORMAT ( ’ EXP( ’ , I , ’ ) = ’ ,F8. 3 , ’ ; ’ ) 

5022 FORMAT ( ’ ETSC AT ( \I F8 . 3 , ’ ; ’ ) 

5023 FORMAT( ’EPSCAT( 1 , 1 F8 . 3 , ' ; ’ ) 

5024 FORMAT ( ’ETHF(' , I , ’ ) = ’ , F8 . 3 , ’ ; ’ ) 

5025 FORMAT (’EPHFC , I , ’ ) = ’ , F8 . 3 , ’ ; ’ ) 

CL0SE(2) 

CL0SE(3) 

CL0SE(4) 

CLOSE (7 ) 

CL0SE(8) 

CL0SE(9) 

CLOSE(IO) 

999 CONTINUE 
STOP 
END 

C********************************************************END OF MAIN PROGRAM. 
C 

C* ********************************************** ********* SUBROUTINE ZMAT . 

C REFERENCE: AN IMPROVED E-FIELD SOLUTION FOR CONDUCTING BOR 
C J.R.MAUTZ AND R.F. HARRINGTON 

C TECHNICAL REPORT TR-80-1 

C ROME AIR DEVELOPMENT CENTER 

C CONTRACT NO. F-30602-79-C-001 1 

SUBROUTINE ZMAT(M1 , M2 , NP , NPHI , NT , RH , ZH , X , A , XT , AT , Z ) 

C 

C COMPUTE THE MM IMPEDANCE MATRIX ELEMENTS. THIS IS FROM MAUTZ AND 
C HARRINGTON (NO CHANGES). 

C 

COMPLEX Z( 100000) ,U1 , U2 ,U3 , U4 , U5 ,U6 , U7 , U8 , U9 , UA , UB , G4A (4) , G5A(4) 
COMPLEX G6A (4 ) , G4B (4) , G5B(4) , G6B (4 ) , H4A , H5A , H6A , H4B , H5B 
COMPLEX H6B,UC,UD,GA(100) ,GB(100) 

DIMENSION RH(400) ,ZH(400) ,X(100) ,A(100) ,XT(4) ,AT(4) 

DIMENSION RS(400) ,ZS(400) ,D(400) ,DR(400 ) , DZ(400) 

DIMENSION DM(400) ,C2(100) ,C3(100) ,R2(4) ,Z2(4) 

DIMENSION C4( 100) ,C5(100) ,C6(100) ,Z7(4) ,R7(4) ,Z8(4) ,R8(4) 

CT=2 . 

CP=. 1 

DO 10 1=2 ,NP 
12=1-1 

RS(I2)=.5*(RH(I) +RH(I2) ) 

ZS(I2)=.5*(ZH(I)+ZH(I2)) 

Dl= . 5* (RH(I)-RH(I2) ) 

D2= . 5* (ZH(I ) -ZH(I2) ) 

D(I2)=SQRT(D1*D1+D2*D2) 

DR(I2)=D1 
DZ(I2) =D2 

IF(RS(I2) .EQ.O. ) RS(I2)=1 . 

DM(I2)=D(I2)/RS(I2) 

10 CONTINUE 



i 0 



M3=M2-M1+1 
M4=M1-1 
PI2=1 .570796 
DO 11 K=1,NPHI 
PH=PI2*(X(K)+1 . ) 
C2(K) =PH*PH 
SN=SIN( ,5*PH) 
C3(K)=4. *SN*SN 
A1=PI2*A(K) 

D4= . 5*A1*C3(K) 

D5=A1*C0S(PH) 

D6=A1*SIN(PH) 

H5=K 

DO 29 M=1 ,M3 

PHM=(M4+M)*PH 

A2=C0S(PHM) 

C4(M5) =D4*A2 
C5(M5)=D5*A2 
C6(M5) =D6*SIN(PHM) 
H5=M5+NPHI 
29 CONTINUE 
11 CONTINUE 
HP=NP-1 
MT=HP-1 
N=MT+MP 
N2N=HT*N 
N2=N*N 
U1=(0 . , .5) 

U2=(0. ,2. ) 

JN=-1-N 

DO 15 JQ=1 , HP 

Kq=2 

iF(jq.Eq.i) Kq=i 
iF(jq.Eq.MP) Kq=3 
Rl=RS(jq) 
zi=zs(jq) 

Di=D(jq) 

D2=DR( jq) 

D3=DZ(jq) 

D4=D2/R1 
D5=DM( jq) 

SV=D2/D1 

CV=D3/D1 

T6=CT*D1 

T62=T6+D1 

T62=T62*T62 

R6=CP*R1 

R62=R6*R6 

DO 12 L= 1 , NT 

R2(L) =R1+D2*XT(L) 

Z2(L)=Z1+D3*XT(L) 
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12 CONTINUE 
U3=D2*U1 
U4=D3*U1 
DO 16 IP=1 , MP 
R3=RS(IP) 

Z3=ZS(IP) 

R4=R1-R3 

Z4=Z1-Z3 

FM=R4*SV+Z4*CV 

PHM=ABS(FM) 

PH=ABS(R4*CV-Z4*SV) 

D6=PH 

IF(PHM.LE.Dl) GO TO 26 
D6=PHM-D1 

D6=SQRT(D6*D6+PH*PH) 

26 IF(IP.EQ.JQ) GO TO 27 
KP=1 

IF(T6.GT.D6) KP=2 
IFCR6.GT.D6) KP=3 
GO TO 28 

27 KP=4 

28 GO TO (41,42,41 ,42) ,KP 
42 DO 40 L=1 ,NT 

D7=R2(L)-R3 
D8=Z2(L)-Z3 
Z7 (L)=D7*D7+D8*D8 
R7 (L)=R3*R2(L) 

Z8(L)= . 25*Z7 (L) 

R8 (L)= . 25*R7 (L) 

40 CONTINUE 

Z4=R4*R4+Z4*Z4 
R4=R3*R1 
RS= . 5*R3*SV 
DO 33 K=1,NPHI 
A1=C3(K) 

RR=Z4+R4*A1 
UA=0 . 

UB=0 . 

IF(RR.LT.T62) GO TO 34 
DO 35 L= 1 , NT 
R=SQRT(Z7 (L)+R7 (L)*A1) 
SN=-SIN(R) 

CS=COS(R) 

UC=AT(L)/R*CMPLX(CS,SN) 
UA=UA+UC 
UB=XT(L)*UC+UB 
35 CONTINUE 
GO TO 36 
34 DO 37 L=1 ,NT 

R=SQRT(Z8(L)+R8(L)*A1) 
SN=-SIN (R) 



CS=COS(R) 

UC=AT(L)/R*SN*CMPLX(-SN ,CS) 

UA=UA+UC 

UB=XT(L)*UC+UB 

37 CONTINUE 
A2=FM+R5*A1 
D9=RR-A2*A2 
R=ABS(A2) 

D7=R-D1 

D8=R+D1 

D6=SQRT(D8*D8+D9) 

R=SQRT(D7*D7+D9) 

IF(D7 . GE. 0 . ) GO TO 38 
Al=ALOG( (D8+D6) * (-D7+R)/D9)/D1 
GO TO 39 

38 Al=ALOG( (D8+D6)/ (D7+R) )/Dl 

39 UA=A1+UA 

UB=A2* (4. / (D6+R)-A1 )/Dl+UB 
36 GA(K)=UA 
GB(K)=UB 
33 CONTINUE 
K1=0 

DO 45 M=1 ,M3 
H4A=0 . 

H5A=0 . 

H6A=0 . 

H4B=0 . 

H5B=0 . 

H6B=0 . 

DO 46 K=1 ,NPHI 

K1=K1+1 

D6=C4(K1) 

D7=C5 (Kl) 

D8=C6(K1) 

UA=GA(K) 

UB=GB(K) 

H4A=D6*UA+H4A 
H5A=D7*UA+H5A 
H6A=D8*UA+H6A 
H4B=D6*UB+H4B 
H5B=D7*UB+H5B 
H6B=D8*UB+H6B 
46 CONTINUE 
G4A(H)=H4A 
G5A(M)=H5A 
G6A(M)=H6A 
G4B(M)=H4B 
G5B(M)=H5B 
G6B(M)=H6B 
45 CONTINUE 

IF(KP . NE.4) GO TO 47 
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A2=D1/(PI2*R1) 

D6=2./D1 
D8=0 . 

DO 63 K=1 ,NPHI 
A1=R4*C2(K) 

R=R4*C3(K) 

IF(R.LT.T62) GO TO 64 
D7=0 . 

DO 65 L= 1 , NT 

D7=D7+AT(L)/SQRT(Z7(L)+A1) 

65 CONTINUE 
GO TO 66 

64 A1=A2/(X(K) + 1. ) 

D7=D6*AL0G(A1+SQRT(1 . +A1*A1) ) 

66 D8=D8+A(K)*D7 
63 CONTINUE 

Al= . 5*A2 
A2=l ./Al 

D8=-PI2*D8+2 ./R1*(BL0G(A2)+A2*BL0G(A1) ) 

DO 67 M=1 ,M3 
G5A(M)=D8+G5A(M) 

67 CONTINUE 
GO TO 47 

41 DO 25 M=1 ,M3 
G4A(M) =0 . 

G5A (M) =0 . 

G6A(M)=0. 

G4B(M)=0. 

G5B(M)=0. 

G6B(M)=0. 

25 CONTINUE 

DO 13 L=1 ,NT 
A1=R2(L) 

R4=A1-R3 

Z4=Z2(L)-Z3 

Z4=R4*R4+Z4*Z4 

R4=R3*A1 

DO 17 K=1,NPHI 

R=SQRT(Z4+R4*C3(K) ) 

SN=-SIN (R) 

CS=COS(R) 

GA(K)=CMPLX(CS ,SN)/R 
17 CONTINUE 
D6=0 . 

IF(R62.LE.Z4) GO TO 51 
DO 62 K=1 ,NPHI 
D6=D6+A(K)/SQRT(Z4+R4*C2(K) ) 

62 CONTINUE 

Z4=3. 141593/SQRT(Z4/R4) 

D6=-PI2*D6+AL0G(Z4+SQRT( 1 . +Z4*Z4) )/SQRT(R4) 
51 A1=AT(L) 
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A2=XT(L) *A1 
K1=0 

DO 30 M= 1 ,H3 
U5=0 . 

U6=0. 

U7 =0 . 

DO 32 K=1,NPHI 
UA=GA(K) 

K1=K1+1 

U5=C4(K1)*UA+U5 
U6=C5(K1)*UA+U6 
U7=C6(K1 ) *UA+U7 
32 CONTINUE 
U6=D6+U6 

G4A(M)=A1*U5+G4A(M) 
G5A(M)=A1*U6+G5A(M) 
G6A(M)=A1*U7+G6A(M) 
G4B(M)=A2*U5+G4B(M) 
G5B(H)=A2*U6+G5B(H) 
G6B(M)=A2*U7+G6B(M) 
30 CONTINUE 
13 CONTINUE 
47 A1=DR(IP) 

UA=A1*U3 

UB=DZ(IP)*U4 

A2=D(IP) 

D6=-A2*D2 

D7=D1*A1 

D8=D1*A2 

JH=JN 

DO 31 H=1 ,M3 

FM=M4+M 

A1=FM*DM(IP) 

H5A=G5A(H) 

H5B=G5B(M) 

H4A=G4A(H)+E5A 

H4B=G4B(M)+H5B 

H6A=G6A(M) 

H6B=G6B(M) 

U7=UA*H5A+UB*H4A 

U8=UA*H5B+UB*H4B 

U5=U7-U8 

U6=U7+U8 

U7=-U1*H4A 

U8=D6*H6A 

U9=D6*H6B-A1*H4A 

UC=D7*(H6A+D4*H6B) 

UD=FM*D5*H4A 

K1=IP+JM 

K2=K1+1 

K3=K1+N 



K4=K2+N 

K5=K2+MT 

K6=K4+MT 

K7=K3+N2N 

K8=K4+N2N 

GO TO (18,20, 19) ,KQ 

18 Z(K6)=U8+U9 
IF(IP.EQ. 1) GO TO 21 
Z(K3)=Z(K3)+U6-U7 
Z(K7)=Z(K7)+UC-UD 
IF(IP.EQ.MP) GO TO 22 

21 Z(K4)=U6+U7 
Z(K8)=UC+UD 
GO TO 22 

19 Z(K5)=Z(K5)+U8-U9 
IF(IP.EQ. 1) GO TO 23 
Z(K1)=Z(K1)+U5+U7 
Z(K7 )=Z(K7 )+UC-UD 
IF(IP.EQ.MP) GO TO 22 

23 Z(K2)=Z(K2)+U5-U7 
Z(K8)=UC+UD 

GO TO 22 

20 Z(K5)=Z(K5)+U8-U9 
Z(K6)=U8+U9 
IF(IP.EQ.l) GO TO 24 
Z(K1)=Z(K1)+US+U7 
Z(K3)=Z(K3) +U6-U7 
Z(K7)=Z(K7)+UC-UD 
IF(IP.Eq.MP) GO TO 22 

24 Z(K2)=Z(K2)+US-U7 
Z(K4) =U6+U7 
Z(K8)=UC+UD 

22 Z(K8+MT )=U2*(D8* (H5A+D4*H5B)-A1*UD) 

JM=JM+N2 

31 CONTINUE 

16 CONTINUE 
JN= JN+N 

15 CONTINUE 
RETURN 
END 

C* *********************************************** ****** **SUBROUTINE SOLVE . 

C REFERENCE : TECHNICAL REPORT TR-80-1 
SUBROUTINE SOLVE(N , IPS ,UL , B , X ) 

C 

C SEE MAUTZ & HARRINGTON FOR DETAILS 

C 

COMPLEX UL(lOOOOO) ,B(800) ,X(800) .SUM 
DIMENSION IPS(800) 

NP=N+ 1 
IP=IPS(1) 

X(l)=B(IP) 
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DO 2 1=2, N 
IP=IPS(I) 

IPB=IP 
IM1=I- 1 
SUM=(0 . ,0. ) 

DO 1 J=1 , IM1 
SUM=SUM+UL(IP) *X( J) 

1 IP=IP+N 

2 X (I ) =B(IPB) -SUM 
K2=N* ( N-l ) 

IP=IPS (N) +K2 
X(N) =X(N)/UL(IP) 

DO 4 IBACK=2 , N 
I=NP-IBACK 
K2=K2-N 
IPI=IPS (I ) +K2 
IP1=I+1 

SUM=(0 . ,0.) 

IP=IPI 

DO 3 J=IP1 , N 
IP=IP+N 

3 SUM=SUM+UL(IP) *X( J ) 

4 X(I)=(X(I)-SUM)/UL(IPI) 

RETURN 

END 

C *** ** ***** *** * *********** ************ ******* ♦♦SUBROUTINE DECOMP 
C REFERENCE: TECHNICAL REPORT TR-80-1 
SUBROUTINE DECOMP (N , IPS, UL) 

C 

C SEE MAUT2 & HARRINGTON FOR DETAILS 
C 

COMPLEX UL( 100000 ), PIVOT, EM 
DIMENSION SCL(800) , IPS (800) 

DO 5 1=1, N 
IPS(I)=I 
RN=0 . 

J1 = I 

DO 2 J=1 , N 

ULM=ABS(REAL(UL( Jl ) ))+ABS(AIMAG(UL(Jl)) ) 

J 1= J 1+N 

IF(RN-ULM) 1,2,2 

1 RN=ULM 

2 CONTINUE 
SCL(I)=1./RN 

5 CONTINUE 
NM1=N-1 
K2=0 

DO 17 K=1 , NM1 
BIG=0 . 

DO 11 I=K,N 
IP=IPS(I) 
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IPK=IP+K2 

SIZE=(ABS(REAL(UL(IPK)))+ABS( AIMAG (UL( IPK ) ) ) ) *SCL( IP) 

IF(SIZE-BIG) 11,11,10 

10 BIG=SIZE 
IPV=I 

11 CONTINUE 
IF(IPV-K) 14,15,14 

14 J=IPS(K) 

IPS(K)=IPS(IPV) 

IPS(IPV)=J 

15 KPP=IPS(K)+K2 
PIVOT=UL(KPP) 

KP1=K+1 

DO 16 I=KP1 , N 
KP=KPP 

IP=IPS ( I ) +K2 
EM=-UL( IP) /PIVOT 
18 UL(IP)=-EM 
DO 16 J=KP1,N 
IP=IP+N 
KP=KP+N 

UL(IP)=UL(IP)+EM*UL(KP) 

16 CONTINUE 
K2=K2+N 

17 CONTINUE 
RETURN 
END 

C ******************************************* ******puNCT ION BLOG 
C REFERENCE: TECHNICAL REPORT TR-80-1 ; (page 56) 

FUNCTION BLOG(X) 

IF(X.GT. . 1) GO TO 1 
X2=X*X 

8L0G=( ( .075*X2- . 1666667) *X2+1 . )*X 
RETURN 

1 8L0G=AL0G(X+SQRT( 1 . +X*X) ) 

RETURN 

END 

C***» **************************************** ********SUBROUTINE PLANE 
C REFERENCE: TECHNICAL REPORT TR-80-1 ; (pages 87-64) 

SUBROUTINE PLANE(M1 , M2 , NP , NT , RH , ZH , XT , AT ,THR , R) 

C 

C PLANE WAVE EXCITATION VECTOR IN THE FAR FIELD. FROM MAUTZ AND HARRINGTON. 
C MODIFIED TO DO ONLY ONE ANGLE AND FREQUENCY PER CALL. 

C ** NEW BESSEL FUNCTION ROUTINE FROM NUM . RECIPES ** 

C 

COMPLEX R( 1600) ,U ,U1 ,UA ,UB , FA ( 1500) ,FB( 1500) , F2A , F2B , FI A , FIB 
COMPLEX U2,U3,U4,U5 

DIMENSION RH(400) ,ZH(400) , XT(4) , AT(4) , R2 (4) ,Z2(4) 

MP=NP-1 

MT=MP-1 

N=MT+MP 
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N2=2*N 

CC=COS(THR) 

SS=SIN (THR) 

U=(0. ,1.) 

Ul=3. 141593*U**M1 

M3=M1+1 

M4=M2+3 

IF(Ml.EQ.O) M3=2 

M5=Ml+2 

M6=M2+2 

DO 12 IP=1 , MP 

K2=IP 

I=IP+1 

DR=.5*(RH(I)-RH(IP)) 

DZ=. 5*(ZH(I)-ZH(IP) ) 
D1=SQRT(DR*DR+DZ*DZ) 
R1=.25*(RH(I)+RH(IP)) 

IF(Rl.EQ.O.) Rl=l . 
Z1=.5*(ZH(I)+ZH(IP)) 

DR= . 5*DR 
D2=DR/R1 
DO 13 L=1 ,NT 
R2(L) =R1+DR*XT(L) 
Z2(L)=Z1+DZ*XT(L) 

13 CONTINUE 
D3=DR*CC 
D4=-DZ*SS 
D5=D1*CC 
DO 23 M=M3 , M4 
FA(M)=0 . 

FB(M)=0 . 

23 CONTINUE 

DO 15 L=1 ,NT 
X=SS*R2(L)*2. 

ARG=Z2(L)*CC 

UA=AT(L)*CMPLX(COS( ARG) ,SIN(ARG)) 
UB=XT(L)*UA 

C THIS LINE REPUCES HARRINGTON'S BLOCK 
DO 25 M=M3 , M4 
BES=BESSJ(M-2,X) 
FA(M)=BES*UA+FA(M) 
FB(M)=BES*UB+FB(M) 

25 CONTINUE 
15 CONTINUE 

IF(Ml.NE.O) GO TO 26 
FA(1)=-FA(3) 

FB(l)=-FB(3) 

26 UA=U1 

DO 27 M=M5,M6 
M7=M-1 
M8=M+ 1 
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F2A=UA*(FA(M8)+FA(M7)) 

F2B=UA*(FB(M8)+FB(M7) ) 

UB=U*UA 

F1A=UB*(FA(M8)-FA(H7)) 

F1B=UB*(FB(M8)-FB(M7) ) 

U4=D4*UA 

U2=D3*F1A+U4*FA(H) 

U3=D3*F1B+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 

R(K2+MT)=-D5*(F2A+D2*F2B) 

R(K5+MT)=D1*(F1A+D2*F1B) 

IF(IP.EQ.l) GO TO 21 
R(K1)=R(K1) +U2-U3 
R(K4) =R(K4) +U4-U5 
IF(IP.EQ.MP) GO TO 22 

21 R(K2)=U2+U3 
R(K5)=U4+U5 

22 K2=K2+N2 
UA=UB 

27 CONTINUE 
12 CONTINUE 
RETURN 
END 

£****************************************************5 UBROUTINE ZLOAD . 
C REFERENCE: COMPUTATION OF RADIATION AND SCATTERING FOR 
C LOADED BODIES OF REVOLUTION 

C HARRINGTON AND MAUTZ 

C REPORT AFCRL-70-0046 

C AIRFORCE CAMBRIDGE RESEARCH LABS 

C CONTRACT NO. F-19628-68-C-0180 

SUBROUTINE ZLOAD (NP , RH , ZH , ZO , Z) 

C 

C COMPUTES IMPEDANCE MATRIX ELEMENTS FOR LOADED BODIES OF REV 
C ZO(I) IS THE SURF IMPEDANCE OF THE ITH SEGMENT (NP-1 SEGMENTS) 

C Z(.) ARE THE IMPEDANCE MATRIX TERMS (TRIDIAGONAL FOR T-T 
C SUBMATRIX; DIAGONAL FOR P-P SUBMATRIX). STORED IN COL VECTOR. 

C 

COMPLEX Cl , C2 , Z0(400) , Z(2400) , X 1 , X2 , X3 , Y1 , Y2 , Y3 , FN(400) 

COMPLEX U1 , U2 ,U3 , XI , YI 

DIMENSION RH(400) ,ZH(400) , RS(400) ,D(400) ,SV(400) 

PI=3. 14159 

MT=NP“2 

MP=NP-1 

N=MT+MP 

DO 10 IP=2 , NP 

II=IP-1 

DR=RH( IP) -RH(II ) 
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DZ=ZH(IP)-ZH(II) 

D(II)=SqRT(DR*DR+DZ*DZ) 

SV (II ) =DR/D ( II ) 

RS(II )= . 5* (RH(IP)+RH(II) ) 

DS=D(II)*SV(II)/2. 

Q1=RS(II)+DS 
Q2=RS(II)-DS 
FN(II)=1 . 

IF((ABS(Q2) .GT. l.E-6) .AND. (ABS(Ql) . GT 1 E-6)) 

* FN( II ) =ALOG (Q1/Q2) 

10 CONTINUE 
L0=MT*3-2 
DO 20 1=1, MP 
C1=PI*Z0(I) 

IF(I.Eq.MP) GO TO 80 
KI=2 

IF(I.Eq.l) KI =1 
IF(I.Eq.MT) KI=3 
11=1+1 

C2=PI*Z0 (II ) 

A=SV (I ) 

IF(ABS(A) .LT. l.E-6) GO TO 41 
Xl=Cl*FN(I)/2. / A 

X2=Cl*2./A*(l.-RS(l)*FN(I)/D(I)/A) 

X3=-X2*RS(I)/D(I)/A 
GO TO 42 

41 CONTINUE 
Xl=Cl/2./RS(I)*D(I) 

X2=(0. ,0.) 

X3=Cl*D(I)/6./RS(I) 

42 CONTINUE 
A=SV (II ) 

IF(ABS(A).LT. l.E-6) GO TO 45 
Yl=C2*FN(II)/2 . / A 

Y2=C2*2./A*(1.-RS(II)*FN(II)/D(II)/A) 
Y3=-Y2*RS(II)/D(II)/A 
GO TO 40 
45 CONTINUE 

Yl=C2/2./RS(II)*D(II) 

Y2=(0 . ,0.) 

Y3=C2*D(II)/6. /RS(II) 

40 CONTINUE 
C 

C DEFINE TRIDIAGONAL ELEMENTS FOR T-T SUBMATRIX (STORED IN COLS) 
C (Ul- DIAG ; U2- LOWER; U3- UPPER) 

C 

XI=X1+X2+X3 

YI=Y1-Y2+Y3 

IF(RH(I).LT. l.E-6) XI=C1/SV(I) 

IF(RH(II) .LT. l.E-6) YI=C2/SV(II) 

U1=XI+YI 
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U2=X1-X3 

U3=Y 1-Y3 

L=2+(I-2)*3 

IF(KI.EQ. 1) L=0 

L1=L+1 

L2=L+2 

L3=L+3 

go to (50,60 ,70) ,ki 
50 Z(L1 )=U1 
Z(L2)=U2 
GO TO 80 
60 Z(L1 )=U3 
Z(L2)=U1 
Z(L3)=U2 
GO TO 80 
70 Z(L1 )=U3 
Z(L2) =U1 

80 Z(L0+I)=2 . *C1*D(I )/RS(I ) 

20 CONTINUE 
RETURN 
END 

SUBROUTINE ZTOT(MT,HP,ZL,Z) 

C 

C ADDS THE SURF IMPEDANCE TERMS TO THE TRIDIAGONAL ELEMENTS OF 
C THE BOR IMPEDANCE MATRIX Z. 

C 

COMPLEX ZL(2400) .Z(IOOOOO) 

N=MT+MP 

M0=MT*3-2 

DO 100 1=1, MP 

L0=MT*N+(I-1)*N+MT 

IF(I.Eq.MP) GO TO 80 

KI=2 

IF(I.EQ.l) KI=1 
IF(I.EQ.MT) KI=3 
L2=(I-1 )*N+I 
L1=L2-1 
L3=L2+1 
M=2+3*(I-2) 

IF(KI.EQ.l) M=0 
M1=M+1 
M2=M+2 
M3=M+3 

go to (50 ,60 , 70) ,ki 
50 Z(L2)=Z(L2)+ZL(M1) 

Z(L3)=Z(L3)+ZL(M2) 

GO TO 80 

60 Z(L1)=Z(L1)+ZL(M1) 

Z(L2)=Z(L2)+ZL(M2) 

Z(L3)=Z(L3)+ZL(M3) 

GO TO 80 
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70 Z(L1)=Z(L1)+ZL(M1) 

Z(L2)=Z(L2) +ZL(H2) 

80 Z(L0+I)=Z(L0+I)+ZL(H0+I) 

100 CONTINUE 
RETURN 
END 

C***** ***************************************** ♦♦SUBROUTINE OGIVE 
C SUBROUTINE : OGIVE 
C DATE : 4 SEPTEMBER 1991 

C PROGRAMMER : R.M. FRANCIS 
C REVISED : 26 JANUARY 1992 

C COMMENTS : 

C THIS SUBROUTINE WILL GENERATE DATA FOR A BODY OF REVOLUTION (BOR) IN 
C THE FORM OF AN OGIVE. 

C DIMENSIONS ARE NORMALIZED TO WAVELENGTH, SEE PAGE 14 FRANCIS THESIS. 

C ZH = Z CO-ORDINATE * 2*PI 
C RH = RADIUS *2*PI 

SUBROUTINE OGIVE(NP , ARAD , ZH , RH , BASE , RS , ZP ) 

C NP = NUMBER OF POINTS ON THE OGIVE SURFACE, MAXIMUM = 1000 
C ZP = ZPRIME , THE POSITION ON Z WHERE THE RADIUS OF CURVATURE STARTS 
C BASE = BASE RADIUS 

C RS = RADIUS OF CURVATURE IN THE RZ,Z PLANE OF THE OGIVE 
INTEGER I , NP , NPBASE 

REAL RH(400) , ZH (400 ) ,ZP, BASE, RS.ZCOORD, RADIUS, AL 
PI=3 . 1415926 

C INPUT THE VARIABLES FOR THE OGIVE, ZP,B,RS,NP 

WRITE(6 , * ) ’ENTER SURFACE CURVATURE (wavelengths)’ 

READ (5 , *) RS 

WRITE(6 , *) ’ENTER ZPRIME, WHERE CURVATURE STARTS (wavelengths)’ 
READ (5 , *) ZP 

WRITE (6 , * ) ’ENTER BASE RADIUS (wavelengths)’ 

READ (5 , *) BASE 
C PERFORM CALCULATIONS 

AL=SQRT(2. *BASE*RS-BASE**2) 

ZMAX= AL + ZP 
ANG=ASIN(AL/RS) 

L=ANINT ( (RS* ANG+ZP) *5) 

NP=2*L+1 



WRITE(6,*) ’NUMBER OF POINTS IS: ’,NP 
DZ= ZMAX/FL0AT(NP-1) 

DO 10 1=1, NP 

ZCOORD=ZMAX- FL0AT(I-1)*DZ 
IF(ZCOORD.EQ.O. )THEN 
ZC00RD=0. 0000000001 
ENDIF 

ZH(I)=2.*PI*ZC00RD 

IF (ZCOORD.LE.ZP) THEN 
RADIUS=BASE 

ELSE 

RADIUS=SqRT(RS**2-(ZC00RD-ZP)**2)+(BASE-RS) 
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IF (RADIUS . EQ . 0 . )THEN 
RADIUS=0. 0000000001 
ENDIF 
ENDIF 

RH ( I ) =2 . *PI*RADIUS 
10 CONTINUE 
RETURN 
END 



(>♦♦♦♦♦♦♦♦♦ ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦SUBROUTINE GENEX 



C SUBROUTINE 
C DATE 
C REVISED 
C PROGRAMMER 



GENEX 

14 JANUARY 1992 
17 February 1992 
R.M. FRANCIS 

SUBROUTINE GENEX (MODE , NP , NT , NPHI , CNPHI , XT , AT , X , A , 

* XP , AP , THS , PHS , ARAD , RHO , ZHO , R) 

C SOURCE ON Z AXIS AT Z=0. THIS VERSION DOES THE ANTENNA INTEGRATION 
C IN RECTANGULAR COORDINATES. CNPHI POINTS USED IN X AND Y. 

C ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

c * 

C t t i * 

C (V) = < W , E > * 

C n l ni * 

C * 

C * 

C phi phi i * 

C (V) = < W , E > * 

C n i ni * 

C * 

♦♦♦♦♦♦♦ ********************* 

C THE EXCITATION VECTOR IS COMPUTED FOR THE GIVEN R, THETA AND PHI 
C COMPONENTS ON SURFACE OF BOR, SPECIFIED IN SUBROUTINE CIRCRTP . 
COMPLEX R(800) , CEXP , PSI , SI , S2 , S3 , S4 , S5 
COMPLEX CIRCR , CIRCT , CIRCP 

DIMENSION RH (400) , ZH (400) , XT(4) , AT(4) ,X(lOO) ,A(100) 

DIMENSION XP(IOO) , AP ( 100) , RH0(400) , ZH0(400) 

INTEGER CNPHI 
PI=3. 141592654 
BK=2 . *PI 



MP=NP- 1 
MT=MP- 1 
N=MT+MP 

C LIMITS ON PHI INTEGRATION 
Pl=(2 . *PI-0 . )/2 . 

P2=(2 . *PI+0 . )/2 . 

DO 10 1=1, NP 
RH (I ) =RHO( I )/BK 
10 ZH(I)=ZHO(I)/BK 
DO 30 IP= 1 , MP 

C QUANTITIES FOR THE FIRST SEGMENT (POSITIVE SLOPE) 
I=IP+1 



II = IP 
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DR=RH(I)-RH(II) 

DZ=ZH(I)-ZH(II) 

Dl=SqRT(DR*DR+DZ*DZ) 

Rl=. 5*(RH(I)+RH(II) ) 

Zl= . 5* (ZH ( I) +ZH( II) ) 

SVP1=DR/D1 

CVP1=DZ/D1 

Vl=ATAN2(DR,DZ+l.E-5) 

C QUANTITIES FOR THE SECOND SEGMENT (NEGATIVE SLOPE) 

C (SKIP IF LAST SEGMENT) 

I=IP+2 

II=IP+1 

DR=RH (I)-RH(II) 

DZ=ZH(I)-ZH(II) 

D2=SQRT(DR*DR+DZ*DZ) 

R2=. 5*(RH(I)+RH(II) ) 

Z2=.5*(ZH(I)+ZH(II)) 

SVP2=DR/D2 

CVP2=DZ/D2 

V2=ATAN2(DR,DZ+l.E-5) 

C BEGIN PHI INTEGRATION: R TERMS IN SI AND S2; THETA TERM IN S3 AND S4; 
C PHI TERM IN S5 . 

S1=(0. ,0. ) 

S2=(0 . ,0. ) 

S3=(0. ,0.) 

S4=(0 . ,0.) 

S5=(0 . ,0. ) 

DO 20 J = 1 ,NPHI . 

PH=P1*X(J)+P2 

PSI=CEXP(CMPLX(0. ,-MODE*PH))*A(J) 

IF(IP.LE.MT) THEN 

C t-CURRENT INTEGRATION FOR THE POSITIVE SLOPE 
C Gauss Quadrature 
DO 13 L=1 ,NT 
TP=Dl*XT(L)/2. 

RHB=R1+TP*SVP1 
ZHB=Z 1+TP*CVP1 
TH=ATAN2(RHB ,ZHB+1 . E-5) 

CC=C0S(V1-TH) 

SS=SIN(V1-TH) 

CALL CIRCRTP(CNPHI ,XP , AP , ARAD ,THS , PHS , 

* PH.RHB ,ZHB ,CIRCR,CIRCT,CIRCP) 

S1=S1+AT(L)*(. 5+TP/Dl)*CC*CIRCR*PSI 
S2=S2+AT(L)*( . 5+TP/Dl)*SS*CIRCT*PSI 
13 CONTINUE 

C t-CURRENT INTEGRATION FOR THE NEGATIVE SLOPE 
C Gauss quadrature 
DO 14 L=1 ,NT 
TP=D2*XT(L)/2. 

RHB=R2+TP*SVP2 

ZHB=Z2+TP*CVP2 
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TH=ATAN2(RHB ,ZHB+1 . E-5) 

SS=SIN ( V2-TH) 

CC=COS( V2-TH) 

CALL CIRCRTP ( CNPHI , XP , AP , ARAD , THS , PHS , 

♦ PH , RHB , ZHB , CIRCR , CIRCT , CIRCP ) 
S3=S3+AT(L)*(. 5-TP/D2) *CC*CIRCR*PSI 
S4=S4+AT(L) * ( . 5-TP/D2) *SS*CIRCT*PSI 

14 CONTINUE 
ENDTF 

C phi-CURRENT INTEGRATION 
DO 15 L=1 , NT 
TP=Dl*XT(L)/2 . 

RHB=R1+TP*SVP1 
ZHB=Z1+TP*CVP1 
TH=ATAN2(RHB , ZHB + 1 .E-5) 

CALL CIRCRTP (CNPHI , XP , AP , ARAD , THS , PHS , 

♦ PH, RHB, ZHB, CIRCR, CIRCT, CIRCP) 
S5=S5+AT(L) *CIRCP*RHB/R1*PSI 

15 CONTINUE 

20 CONTINUE 

C COMPONENTS ARE STORED IN A COLUMN VECTOR: VT(1,MT), VP(MT+1,N) 
IF(IP.LE.MT) THEN 

R(IP) =(Sl+S2)*Dl/2 . *Pl+(S3+S4)*D2/2 . *P1 
ENDIF 

R(IP+MT)=S5*Pl*Dl/2 . 

30 CONTINUE 

RETURN 
END 



C* *♦♦****♦*♦*♦*******♦******♦♦♦*♦**♦♦********** ♦♦♦♦SUBROUTINE CIRCRTP 



SUBROUTINE 

DATE 

REVISED 

PROGRAMMER 

2ND REVISION 

PROGRAMMER 



CIRCRTP 
1 OCTOBER 1991 

26 February 1992 
R.M. FRANCIS 

27 October 1992 
K. A. KLOPP 



COMMENTS : THIS SUBROUTINE WILL RIGOROUSLY CALCULATE THE ELECTRIC FIELD 

FOR A CIRCULAR APERTURE. THE FIELD IS CALCULATED AT THE COORDINATES 
SPECIFIED BY RH(.) AND ZH( . ) THE APERTURE IS LOCATED AT Z = 0, 

AND SCANNED TO A DIRECTION (THS, PHS). THE SUBROUTINE OGIVE IS THE 
SOURCE OF THE GEOMETRIC DATA REQUIRED BY CIRC TO PERFORM COMPUTATIONS. 
ALL PHYSICAL DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 



C +ant +ant * 
C * 
Ci \ \ xl * 

C E = -jn/(4*pi*k) /_ /_ e ( )*TAPER(rp) * 

C x=0 y=0 * 
C * 
c***************************************************************************** 
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SUBROUTINE CIRCRTP(CNPHI , XP , AP , ARAD , THS , PHS , 

* PHI , RZ,Z,CIRCR,CIRCT,CIRCP) 

INTEGER CNPHI 
DIMENSION XP(IOO) ,AP(100) 

REAL K 

COMPLEX J, JK,G1,G2,X1,C0N,CC,SUMX,SUMY,SUMZ 
COMPLEX CIRCR ,CIRCT ,CIRCP 
PI=3 . 141592654 
ETA=120.*PI 

C SET FIELD COMPONENTS TO ZERO IN THE REAR HEMISPHERE 
CIRCR=(0. ,0. ) 

CIRCT=(0. ,0. ) 

CIRCP=(0 . ,0. ) 

C Z OF THE OBSERVATION POINT IS >= 0 IN THE FORWARD HEMISPHERE 
IFCZ.GE.O.) THEN 
SUMX=(0. 0,0.0) 

SUMY=(0. 0,0.0) 

SUMZ=(0. 0,0.0) 

K=6. 283185307 
J = (0.0, 1.0) 

JK= (0.0,6.283185307) 

C0N=-J/(4.*PI) 

C OMIT THE FACTOR OF ETA SINCE IT IS NOT INCLUDED IN SUBROUTINE ZMAT . 
C ARAD**2 IS FOR SCALING OF GAUSSIAN INTEGRATION. READJUST SCALING 
C TO AGREE WITH OLD SUBROUTINE SPHERE (ARBITRARY) 

CC= CON 

C OUTER LOOP: INTEGRATE OVER X, -ARAD < X < ARAD. 

C INNER LOOP: INTEGRATE OVER Y, -ARAD < Y < ARAD. 

C EVALUATE ANTENNA FIELD AT RZ,Z: 

XO=RZ*COS(PHI) 

Y0=RZ*SIN(PHI) 

STS=SIN(-THS) 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 

C INTEGRATE IN X,Y COORDINATES 
DO 50 M=l, CNPHI 
X=ARAD*XP(M) 

DO 60 N=l, CNPHI 
Y=ARAD*XP(N) 

RP=SQRT(X**2+Y**2) 

C SKIP INTEGRATION POINTS BEYOND THE ANTENNA RADIUS 
IF(RP.LE.ARAD) THEN 
PHIPRIME=ATAN2(Y ,X+1 . E-10) 

R=SQRT( ((RZ-RP)**2+Z**2)+(4*RZ*RP*SIN( (PHI-PHI PRIME) /2) **2) ) 
G1=(((K*R)**2)-1-(JK*R))/R**3 
G2=(3+(3*JK*R)-(K*R)**2)/R**5 
X1=JK*(RP*C0S(PHIPRIME-PHS)*STS-R) 

Y1=X0-X 
X 2=Y1**2 
Y2=Y0-Y 

SUMX=SUMX+(CC*AP(M)*AP(N)*(G1+(X2*G2))*CEXP(X1))*TAPER(RP) 
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SUMY=SUMY+(CC*AP(M) *AP(N)*Y1*Y2*G2*CEXP(X1) )*TAPER(RP) 
SUMZ=SUMZ+(CC*AP(M)*AP(N)*Y1*Z*G2*CEXP(X1))*TAPER(RP) 
ENDIF 

60 CONTINUE 

50 CONTINUE 

C ASSUME AN ELEMENT FACTOR, ELMT 
c ELMT=SQRT ( ABS ( C ) ) 

ELMT=1 . 

CIRCT= (C*COS( PHI )*SUMX+C*SIN( PHI) *SUMY-S*SUMZ)* ELMT 

CIRCR= (S*COS (PHI ) *SUMX+S*SIN ( PHI) *SUMY+C*SUMZ) *ELMT 

CIRCP= ( -SIN (PHI ) *SUMX+COS (PHI ) *SUMY ) *ELMT 

ENDIF 

RETURN 

END 



C*** *************************************** ****SUBROUTINE TESTSPHERE 
C SUBROUTINE : TESTSPHERE 

C DATE 21 FEBRUARY 1992 

C PROGRAMMER R. M. FRANCIS 

C COMMENTS GENERATES A SPHERE, WITH PLOTTED POINTS 

C AT EQUAL INTERVALS IN THETA. 

SUBROUTINE TESTSPHERE ( NP , ZH , RH , BASE , RS , ZP ) 

INTEGER NP , I 

REAL ZH (400) ,RH(400) , BASE , RS , ZP , SPRAD 
PI=3. 1415926 
BK=2 . *PI 



ZP=0.0 

WRITE(6,*) ‘ENTER NUMBER OF POINTS (NP):’ 

READ (5 , *)NP 

WRITE(6,*) ‘ENTER SPHERE RADIUS’ 

RE AD (5,*) SPRAD 
BASE=SPRAD 
RS=SPRAD 
DO 1241 1=1, NP 

ANGLE=PI*FL0AT(I-1 )/(2. *FL0AT(NP-1 ) ) 

ZH(I)=BK*SPRAD*COS (ANGLE) 

RH(I )=BK*SPRAD*SIN( ANGLE) 

IF(RH(I).EQ.O.)THEN 

RH(I)=0. 0000000001 
ENDIF 

1241 CONTINUE 
RETURN 
END 

C***«* ********************************* ************SUBROUTINE CONE 
C GENERATES THE BOR GEOMETRY FOR A CONE. 

SUBROUTINE CONE(NP , ZH , RH , BASE , RS , ZP) 

REAL ZH(400) ,RH(400) ,BASE,RS,ZP 
REAL HA,ANG,ZMAX,DZH 
INTEGER NP 
PI=3. 1415926 



BK=2 . *PI 

WRITE (6,*) ‘ENTER THE CONE HALF ANGLE IN DEGREES’ 
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READ (5 , *)HA 

IF(HA . GT. 90 . )THEN 

WRITE (6,*) ’RANGE OF ANGLE: 0< ANGLE<90 ’ 

WRITE(6,*) ’RANGE EXCEEDED’ 

GOTO 10 
ENDIF 

ANG=HA*PI/180. 

WRITE(6,*) ’ENTER THE MAXIMUM Z’ 

READ(5 , *)ZMAX 

WRITE(6,*) ’ENTER THE NUMBER OF POINTS (NP)’ 

READ (5 , *)NP 

ZH(l)=ZMAX*BK 

RH(1)=0. 

DZH=ZH( 1)/FLQAT(NP-1 . ) 

DO 10 1=1, NP 

IF(ZH(I).EQ.O. )THEN 
ZH(I)=. 00000001 
ENDIF 

ZH(I+1)=ZH(I)-DZH 
RH(I+1)=RH(I)+(DZH*SIN(ANG) ) 

10 CONTINUE 

BASE=RH(NP)/BK 

ZP=0. 

RS=0 . 

RETURN 

END 

C**********************************************SUBROUTINE DISK 
SUBROUTINE DISK (NP , ZH , RH , BASE , RS , ZP ) 

REAL ZH(400) ,RH(400) , BASE , RS , ZP , RMAX 
INTEGER NP 

WRITE(6,*) ’ENTER THE RADIUS OF THE DISK’ 

READ(5 , *)RMAX 

WRITE (6,*) ’ENTER THE DISTANCE TO THE DISK’ 

READ (5 , *)ZP 

WRITE (6,*) ’ENTER THE NUMBER OF POINTS (NP)’ 

READ(5,*)NP 
PI=3. 1415926 
BK=2 . *PI 

DRH=RMAX*BK/FL0AT(NP-1 . ) 

DO 10 1=1 ,NP 
ZH(I)=ZP*BK 
RH(I)=DRH*FL0AT(I-1 . ) 

10 CONTINUE 
BASE=RMAX 
RS=0 . 

RETURN 

END 

C* ************************************ *********SUBROUTINE PARAB 
SUBROUTINE P ARAB (NP , ZH ,RH , DM , FOD ,ZP ) 

DIMENSION RH(400) ,ZH(400) 

WRITE (6 , *) ’ENTER DIAMETER AND F/D RATIO’ 
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READ ( 5 , * ) DM , FOD 

WRITE (6 , * ) ’FEED DISTANCE FROM FOCUS (+ IS NEARER)’ 
READ (5 , *) ZP 

WRITE (6, «) ’ENTER NUMBER OF GENERATING POINTS’ 

READ (5 , *) NP 
BK=2.*3. 14159 
FM=FOD*DM 

PHIV=2 . *ATAN( 1 . /4 . /FOD) 

DO 24 1=1, NP 

TH=FLO AT( I- 1 ) *PHIV/FLOAT(NP- 1 ) 

RM=2. *FM/(1 ,+COS(TH))*BK 

ZH(l)=RM*COS(TH)+ZP*BK 

RH(I)=RM*SIN(TH) 

24 CONTINUE 
RETURN 
END 

C**** ******* ********* *************** ******* FUNCTION TAPER 
FUNCTION TAPER(RHO) 

C SPECIFY ANTENNA ILLUMINATION FUNCTION. REAL FUNCTION OF 
C RHO ONLY 

TAPER=1 . 

RETURN 

END 

C*** ***************************** *********SUBROUTINE ANTFF 
SUBROUTINE ANTFF (THETA , PHI , THETAS , PHIS , ARAD , ETF , EPF) 

C COMPUTES CLOSED FORM ANTENNA PATTERN IN THE FAR FIELD. 

C ETF AND EPF ARE RETURNED. ANGLES ARE IN RADIANS. 

C PRESENT ENTRY IS FOR: 

C 

C UNIFORM ILLUMINATION SCANNED TO A DIRECTION (THETAS , PHIS) 
C 

COMPLEX ETF , EPF , JK , SCL , EB 
JK=(0. 0,6. 283185307) 

BK=6. 283185307 
PI=BK/2 . 

C SET RADIATION PATTERN TO ZERO IN THE REAR HEMISPHERE 
IF (COS (THETA) . LT.O. ) THEN 
EB=(0. ,0. ) 

ELSE 

BB=S IN ( THETA ) *SIN( PHI) -SIN (THETAS )*SIN( PHIS) 

AA=SIN (THETA )*COS( PHI )-SIN(THETAS)*COS (PHIS) 
ARG=BK*ARAD*SQRT(AA**2+BB**2) 

SCL=-(0. , 1 . ) *2 . *PI**2 
IF(ABS ( ARG) ,LT.1.E-5)THEN 
EB=SCL/2. 

ELSE 

EB=SCL*(BESSJ1( ARG)/ARG) 

ENDIF 

ENDIF 

C EB IS THE X COMPONENT; GET ETHETA (ETF) AND EPHI (EPF) 

C ASSUME AN ELEMENT FACTOR, ELMT 
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ELMT=SQRT( ABS(C0S (THETA ) ) ) 

ELMT=1 . 

ETF=EB*COS (THETA) * COS (PHI ) *ELMT 
EPF=-EB*SIN(PHI ) *ELMT 
RETURN 
END 

FUNCTION BESS J(N , X) 

C RETURNS THE BESSEL FUNCTION B OF ORDER N (>1) AND REAL 
C ARGUMENT X. 

PARAMETER ( IACC=40 , BIGNO=l . ElO , BIGNI= 1 . E- 10) 
IF(N.EQ.O) THEN 
C IF N=0 CALL BESS JO 

BESS J=BESS JO (X) 

ELSE IF(N.EQ.l) THEN 
C IF N=1 CALL BESS J 1 

BESS J= BESS J 1 (X ) 

ELSE 

C IF N>1 USE RECURSION 
BESS J=0 . 

IF(X . NE . 0 . ) THEN 
T0X=2./X 

IF(X . GT.FLOAT(N) ) THEN 
B JM=BESSJO (X) 

B J=BESS J1 (X) 

DO 11 J=1 ,N-1 
B JP= J*TOX*B J-B JM 
B JM=BJ 
B J=B JP 

11 CONTINUE 
BESS J=B J 
ELSE 

M=2*( ( N+INT(SQRT(FLOAT( IACC*N) ) ) )/2) 

BESS J=0 . 

JSUM=0 . 

SUM=0 . 

B JP=0 . 

B J=1 . 

DO 12 J=M , 1 , -1 
B JM= J*TOX*B J-B JP 
B JP=B J 
B J=B JM 

IF(ABS(BJ) .GT.BIGNO) THEN 

B J=B J*BIGNI 

B JP=B JP*BIGNI 

BESS J=BESS J*BIGNI 

SUM=SUM*BIGNI 

ENDIF 

IF(JSUM.NE.O) SUM=SUM+BJ 
JSUM= 1- JSUM 
IF(J.EQ.N) BESS J=BJP 

12 CONTINUE 



96 



SUM=2 . *SUM-BJ 

BESS J=BESSJ /SUM 

ENDIF 

ENDIF 

ENDIF 

RETURN 

END 

FUNCTION BESSJO(X) 

C 

C BESSEL FUNCTION OF 0 ORDER, REAL ARGUMENT X 
C (SEE ’NUMERICAL RECIPES’, P.172) 

C 

REAL*8 Y,P1 ,P2,P3,P4,P5,Ql,Q2,q3,Q4,Q5,Rl ,R2,R3,R4,R5,R6, 

* S1,S2,S3,S4,S5,S6 

DATA PI , P2 , P3 , P4 , P5/1 . DO , - . 109862827D-2 , . 27345 10407D-4 , 

* -.2073370639D-5, . 2093887211D-6/ 

DATA qi,q2,q3,Q4,q5/-. 1562499995D-1, . 1430488765D-3 , 

* -.6911 147651D-5 , . 7621095161D-6 . 934945152D-7/ 

DATA R1 , R2 ,R3 ,R4,R5, R6/57568490574 . DO , - 13362590354 . DO , 

* 651619640 . 7D0 ,-l 1214424 . 18D0 , 77392 . 33017D0 , -184 . 9052456D0/ 
DATA S 1 , S2 , S3 , S4 , S5 , S6/5756849041 1 . DO , 1029532985 . DO , 

* 9494680. 7 18D0, 59272. 64853D0, 267. 8532712D0, 1 .DO/ 

IFCX.EQ.O.) THEN 

BESSJ0=1 . 

ELSE IF(ABS(X) .LT.8. ) THEN 
Y=X**2 

BESS J0= (Rl+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6) ) ) ) )/ 

* (S1+Y*(S2+Y*CS3+Y*(S4+Y*(S5+Y*S6) ) ) ) ) 

ELSE 

AX=ABS(X) 

Z=8 . / AX 
Y=Z**2 

XX=AX-. 785398164 

BESS JO=SqRT(. 636619772/AX) *(COS (XX)* (PI +Y*(P2+Y* (P3+ 

* Y*(P4+Y*P5) ) ) )-Z*SIN(XX)*(qi+Y*(q2+Y*(q3+ 

* Y*(q4+Y*q5))))) 

ENDIF 

RETURN 

END 

FUNCTION BESSJl(X) 

C 

C BESSEL FUNCTION B OF ORDER 1, REAL ARGUMENT X 
C (SEE ’NUMERICAL RECIPES’, P.173) 

C 

REAL*8 Y,Pl,P2,P3,P4,P5,qi,q2,q3,q4,q5,Rl,R2,R3,R4,R5,R6, 

* SI , S2 , S3 , S4 ,S5 , S6 

DATA PI , P2 , P3 ,P4 , P5/1 . DO , . 183 105D-2 , - . 3516396496D-4 , 

* . 2457520 174D-5 , - . 203370 19D-6/ 

DATA q 1 , q2 , q3 , q4 , q5/ . 04687499995D0 , - . 2002690873D-3 , 

* . 8449199096D-5 , - . 99228987D-6 , . 105787412D-6/ 

DATA R1 , R2 , R3 , R4 , R5 , R6/723626 14232 . DO , -7895059235 . DO , 
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* 242396853.100,-2972611.43900,15704.482600,-30.1603660600/ 
DATA SI , S2 , S3 , S4 , S5 , S6/144725228442 . DO , 2300535178 . DO , 

* 18583304 . 74D0 , 99447 . 4339400 , 376 . 9991397D0 , 1 . DO/ 
IFCX.Eq.O.) THEN 

BESS J 1=0 . 

ELSE IF(ABS(X) .LT.8. ) THEN 
Y=X**2 

BESS J 1=X* (Rl+Y* (R2+Y*(R3+Y*(R4+Y* (R5+Y*R6) ) ) ) )/ 

* (S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) ) ) ) ) 

ELSE 

AX=ABS(X) 

Z=8 . / AX 
Y=Z**2 

XX=AX-2. 356194491 

BESSJ1=SQRT( .63661 9772/AX)* (COS (XX) *(P1+Y*(P2+Y*(P3+ 

* Y*(P4+Y*P5))))-Z*SIN(XX)*(qi+Y*(q2+Y*(q3+ 

* Y*(q4+Y*q5)))))*SIGN(l. ,X) 

ENDIF 

RETURN 

END 
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APPENDIX C. 



GAIN.F LISTING 



Material on the following pages is a listing of the program GAIN.F described 
in chapter III. 
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c******** *************************************************************** 

c 

C PROGRAM : gain.f 

C DATE : 17 November 1992 

C PROGRAMMERS : Dr. D. Jenn, LT K. Klopp 
C 

C Computes the gain of a circular aperature antenna radiating through 
C a radome consisting of a body of revolution, by numerically 
C integrating electric field surrounding system. E-Field is generated by 
C currents distributed on a surface of the body of revolution. The 
C currents, and their locations were previously determined using the 
C program "radome.f”, and stored in the file M f 77curcoef sA or B." 

C 

C INPUT < f 77 curcoef sA/B 
C 

Qt**************************************************************v******* 

CHARACTER ITH 

CHARACTER* 12 CC 

CHARACTER* 10 GO 

CHARACTER*8 GMT , GNPHI , CGP 

CHARACTER* 14 TPTS.PPTS.PHIPTS 

COMPLEX R(1600) .LEADTERM , ETT , EPT 

COMPLEX EP , ET , ETF , EPF , ETCOMB , EPCOMB 

COMPLEX EXP1 , EXP2 , CON JG , CEXP , CMPLX , CT (30000 ) , IMP , JK 

COMPLEX U , UC , ET 1 , EP 1 , ET2 , EP2 , ZLO (400) 

C DIMENSION RH(400) ,ZH(400) , IPS (800) 

DIMENSION RH(400) ,ZH(400) 

DIMENSION XT(4) ,AT(4) ,A(100) ,X(100) ,XP(100) , AP(100) 

DIMENSION Q(30) , S(30) 

DIMENSION THRR(1000, 1000) ,PHRR( 1000) 

C DIMENSION ETSCAT(500) ,EPSCAT(500) 

C DIMENSION EXP(SOO) ,ANG(500) ,ECP(500) 

INTEGER CNPHI, SELECTION 
DATA PI/3.1415926/ 

DATA RR/ 1000/ 

Rad=PI/ 180 . 

BK=2 . *PI 
U=(0. ,1.) 

U0=(0. ,0.) 

UC=-U/4./PI 

JK=(0. 0,6. 283185307) 

C*****ENTER A LETTER TO INDICATE WHICH ITERATION THIS IS***************** 
C***** (ALL DATA FILES ARE APPENDED WITH THIS LETTER) *********** 

C***** (.m FILES ARE FOR GENERATING PLOTS IN MATLAB) *********** 

WRITE(6,*) ’ENTER A LETTER TO INDICATE WHICH ITERATION THIS IS’ 

READ(5 , *)ITH 

CC= ’ curcoef sdat ’ //ITH 

G0= ’gainout ’//ITH// ’ .m’ 

C************************************************* #****#♦*#*******#****** 

C 

C READ in the following data from file called f77curcoefs: 
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c 

C iteration, ITH 

C geometry of BOR, SELECTION 

C radius of base of BOR, BASE 

C total points, NP 

C number of t-intervals, MT 

C phi-intervals , MP 

C sum of t and phi intervals, N 

C number of rows of currents, NROW 

C total number of modes , MODES 

C number of sub vectors, NBLOCK 

C antenna scan angle in theta, THS 

C antenna scan angle in phi, PHS 

C antenna radius, ARAD 

C observation angle, PHIO 

C complex radome impedance, IMP 

C zpnme, ZP 

C surface radius, RS 

C closed form antenna pattern, IFF=0 

C 



C files where Gauss integration weights and abscissas are stored 

C 

C coordinates of surface points *2PI/lamda, ZH(I), and RH(I) 

C 

£**************************★****★★****★***************★★**★*★************* 

0PEN(1 ,f ile=CC) 

WRITE(6,*) ’FILE OPENED IS ’ ,CC 
READ ( 1 , * ) ITH , SELECTION , BASE , NP , MT , MP , N 
READ(1 ,*) NROW.MODES.NBLOCK.THS.PHS.ARAD.PHIO 
READ(1,*) IMP,ZP,RS,IFF 
c write (6, *) IMP , ZP ,RS , IFF 

READ(1,*) (RH(L) ,L=1 ,NP) 
c write (6,*) (RH(L) ,L=1,NP) 

READ(1,*) (ZH(L) ,L=1 ,NP) 
c write(6,*) (ZH(L) ,L=1 ,NP) 

READ(1,*) (ZLO(L) ,L=1,NP) 

C write (6,*) (ZLO(L) ,L=1,NP) 

£************************************************************************** 

c 

C READ in current coefficients, CT(I) from file curcoef sdatA/B 
C (This is one one collumn vector.) 

C 

C************************************************************************** 

READ(1,*) (CT(L) ,L=1 ,NROW) 

C write (6,*) (CT(L) ,L=1 ,NR0W) 

READ ( 1 , *) GNT,GNPHI,CGP 

WRITEC6,*) ’RUN CODE READ FROM DISC : ’,ITH 

THETAS =THS* RAD 

PHIS=PHS*RAD 

MHI=MQDES+1 

CLOSE ( 1 ) 
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c*************************************************«************»»*»**»**»****** 

c 

C Read in weights A(I), and Abscissas X(I) for the following Gauss- 
C quadrature Integrations: 

C file wts absc 

C 

C BOR integration in t direction: TPTS = gaus/gaus#; XT(I),AT(I) 

C BOR integration in phi direction: NPHI = gaus/gaus#; X (I), A (I) 

C Antenna integration in x,y direction: PHIPTS = gaus/gaus#; XP(I),AP(I) 

C 

Cm ******************* *************** ******************************************* 

TPTS= ’ gaus/ ’ / / GNT 

PPTS= ’gaus/ ’//GNPHI 

PHIPTS= ’gaus/ ’ / /CGP 

OPEN ( 2 , FILE=TPTS , STATUS= ’ OLD ’ ) 

OPEN ( 3 , FILE=PPTS , STATUS= ’ OLD ’ ) 

OPEN (8 , FILE=PHIPTS , STATUS^ ’OLD ' ) 

READ (2,*) NT 
WRITE(6,*) ’NT = ’ ,NT 
IF(NT.GT.4)THEN 

WRITE(6,*) ’MAXIMUM NUMBER OF POINTS(NT) IS 4’ 

GOTO 999 
ENDIF 

READ(3,*)NPHI 
WRITE(6,*) ’NPHI = ’ ,NPHI 
IF (NPHI .GT . 200) THEN 

WRITE(6,*) ’MAXIMUM NUMBER OF POINTS(NPHI) IS 200’ 

GOTO 999 
ENDIF 

READ (8 , *)CNPHI 
WRITE(6,*) ’CNPHI = ’, CNPHI 
IF (CNPHI . GT . 100) THEN 

WRITE(6,*) ’MAXIMUM NUMBER OF POINTS(CNPHI ) IS 100’ 

GOTO 999 
ENDIF 

C LOAD THE WEIGHTS AND ABSCISSAS IN THE VECTORS. 

DO 1 M= 1 , NT 

READ (2 , * ,END=1) XT(M) , AT(M) 

WRITE(6 , * )XT(M) , AT(M) 

1 CONTINUE 

DO 2 M=1 .NPHI 

READ(3,*,END=2)X(M) ,A(M) 

WRITE(6 , *) X(M) , A (M) 

2 CONTINUE 

DO 4 M=l, CNPHI 

READ ( 8 , * , END=4 ) XP ( M ) ,AP(M) 

WRITE(6,*)XP(M) , AP(M) 

4 CONTINUE 
CL0SE(2) 

CLOSE (3) 

CL0SE(8) 
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c****************************************************************************** 

c 

C Create an output file called ’ gainoutA/B .m ’ , and output data on geometry etc. 
C 

C****************************************************************************** 

OPEN ( 7 , FILE=G0 ) 

WRITE(7 , 2000 ) * GEOMETRY = * , SELECTION , ’ NP = \NP 
ft, * MT = 1 ,MT, 1 MP = * ,MP, * N = ’ ,N 

ft, 1 MODES = \ MODES,’ NBLOCK = 1 , NBLOCK 

ft, ’ PHI = * ,PHI 

VRITE(7 , 200 1 ) * BASE = * , BASE , * ARAD = * , ARAD , * THETASCAN = ' ,THS 
ft, * PHISCAN = *,PHS 

WRITE(7 , *) J NPHI = \NPHI 

2000 FORMAT (//5X, RADOME GAIN CALCULATION *///. ’ // , 10 (/ »*/. \A13,I4)) 

2001 FORMAT ( 3 ( / * 7, \ A13 ,F9 . 2) , /) 

c 

C Gain Computation: 

C 4 pi * Max Radiation Intensity 

C G = 

C Integration Over Closed Surface ( Radiation 

C Intensity * Beam Solid Angle ) 

C 

C Radiation Intensity = ( 1/ (2*eta) ) * I E(t,phi,r) 1**2 
C 

C Since E is in the far field, need to only determine Etheta, Ephi . 

C 

C | El * *2 = |E scattered I **2 + |E antenna|**2 
C 

C |E scatteredl = Etheta + Ephi 
C 

C******^*****^************:*^*************************************************** 

C**Enter # divisions for integration in theta, and phi 

C directions. (This allows easy variation of the number of integration steps 
C without having to recompute a different set of weights and abscissas for 
C Gauss Integration.) 

WRITE(6,*) ’ENTER NUMBER OF DIVISIONS IN THETA DIRECTION', 

& ’ NDIVT = ?’ 

READ(S , * )NDIVT 

WRITE(6,*) 'NDIVT = NDIVT 

WRITE (7 , *) NDIVT = ’.NDIVT 

WRITE (6,*) 'ENTER NUMBER OF DIVISIONS IN PHI DIRECTION, NDIVP = ?’ 

READ(5 , * ) NDIVP 
WRITE(6,*) 'NDIVP = ’.NDIVP 
WRITE(7,*) ’*/. NDIVP = ', NDIVP 
kk=NDIVT 

C do 500 NDIVT=1 , kk 

NDIVP=NDIVT 

C***********************************integration intervals 0<theta<180, 0<phi<360 
THETA 1=0 
THETA2=180 
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S1=THETA1 

S2=THETA2 

Sl=Sl*Rad 

S2=S2*Rad 

PHI 1=0 

PHI2=360 

Q1=PHI 1 

Q2=PHI2 

Ql=Ql*Rad 

Q2=Q2*Rad 

C***********************************theta integration angles 
DS=(S2-S1) /FLOAT (NDIVT) 

DO 200 1=1 , NDIVT +1 
S(I)=(I-l)*DS 

200 CONTINUE 

C***********************************phi integration angles 
DQ=(Q2-Q1)/FL0AT(NDIVP) 

DO 201 1=1 , NDIVP+l 
q(I)=(I-l)*DQ 

201 CONTINUE 
SUM=0 . 

SUMN0NE=0 . 

UMAX=0 . 

UMAXN0NE=0 . 

EMAX=0 . 

EMAXN0NE=0 . 

C************************v phi integration of rad intensity * beam solid angle V 
DO 300 J J=1 , NDIVP 
Pl=DQ/2. 

P2=(Q(JJ+l)+Q(JJ))/2. 

0************+********************V phi integration determining Escattered V 
C V within each interval of phi, NDIVP V 

DO 300 J=1 , NPHI 
PHR=P1*X( J )+P2 

0****************v theta integration of rad intensity * beam solid angle V 
DO 301 11=1 , NDIVT 
Tl=DS/2 . 

T2=(S(II+1 )+S(II))/2. 

C**************************V theta integration determining Escattered V 
C V within each interval of theta, NDIVT V 

DO 301 1=1, NPHI 
THR=T1*X(I ) +T2 
IF(THR.LT. PI) GO TO 302 
THR=(BK-THR) 

PHR=PHR+PI 

302 CONTINUE 

(>**************** **********************************£ e am Solid Angle 
STRA=SIN(THR)*A(I) *A( J) 

KJ=( J J-l) *NPHI+J 
KI=(II-1)*NPHI+I 
THRR(KJ ,KI)=THR*180/PI 
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PHRR(KJ)=PHR*180/PI 
SNP=SIN (PHR) 

CNP=COS (PHR) 

*********************************** ******************************* ******* 



c 


_ _ 






_ 






c 


1 1 






1 




1 1 


c 


1 s I 






1 _ t, theta 


_ phi, theta 


1 _ t 1 


c 


1 E I 


-jkr 


MHI 


1 R 


R 


1 I 1 


c 


1 t I 


- jhe 


— 


1 n 


n 


1 n | jnphi 


c 


1 1 = 




\ 


1 




1 1 e 


c 


1 s 1 


4pir 


/„ 


1 _ t , phi 


_ phi, phi 


1 _ phil 


c 


1 E I 




n= 1 


1 R 


R 


1 I 1 


c 


1 phil 






1 n 


n 


1 n 1 


c 


1 _ _ 1 






1 _ 







CV********************************see page 25. H,E, and Combined Field********V 
C -jkr 

C -jhe 

C Leading Term : eta (h) is cancelled by its inverse in CT (I), 

C 4pir r dependence is ignored, and the whole Terra is 

C scaled by 2pi for gauss quadrature integration 

C 

C Leading Terra is reduced to (-j/4pi)/2pi 

C*****************************************************************************y 

C LEADTERM=UC/BK 

LEADTERM=UC 
ET=(0 , 0) 

EP=(0 ,0) 

ET1=(0 . ,0 . ) 

EP1 = (0 . ,0. ) 

ET2= (0 . ,0. ) 

EP2= (0 . ,0. ) 

C**************************************** ********************************* ****y 

C MHI 

C 

C \ 

C /„ 

C n=l 

c*****************************************************************************v 

DO 305 M=1 , MHI 
NM=M-1 

C*********************************** t*******************^****** jnphi 
C e 

EXP 1 =CEXP ( CMPLX ( 0 . , FLO AT ( NM ) *PHR) ) 

EXP2=C0NJG(EXP1 ) 

C*************determine plane wave excitation vector in far field, R 
CALL PLANE(HM,HM , NP , NT, RH ,ZH ,XT , AT,THR,R) 

NT0P1=H0DES-KH 
NT0P2=NBL0CK- ( NTOP 1 + 1 ) 

NS2=NT0P1*N 

NS1=NT0P2*N 

C*** ************************************************ **************** 
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c I _ t , theta I I _ t I jnphi I _ t,phi I I _ t I jnphi 

C I R I . I I I e , and I R I . | I |e 

C In I I n I In I I n I 

c******************************************************************* 

DO 306 L=1 , MT 

C***********************************************posit ive mode +n 
ET1=ET1+R(L) *CT(L+NS1 )*EXP1 
EP1=EP1+R(L+N)*CT(L+NS1)*EXP1 

C* ************* ********** ******************** ***negative mode -n 
IF(NM.Eq.O) GO TO 306 
ET1=ET1+R(L)*CT(L+NS2)*EXP2 
EP1=EP1-R(L+N)*CT(L+NS2)*EXP2 

306 CONTINUE 

C****** ************************************************************* 

C I _ phi, theta I I _ t I jnphi I _ phi, phi I I _ t I jnphi 

C|R I . I I |e , and I R I . I I |e 

Cln I I n I In llnl 

C****** ******************* ****************************************** 

DO 307 L=1 , HP 

C***** ************************************** ****posit ive mode +n 
ET2=ET2+R(L+MT)*CT(L+NS1+MT)*EXP1 
EP2=EP2+R(L+MT+N)*CT(L+NS1+MT )*EXP1 
C***********************************************negative mode -n 
IF(NM.EQ.O) GO TO 307 
ET2=ET2-R(L+MT)*CT (L+NS2+MT) *EXP2 
EP2=EP2+R(L+MT+N)*CT (L+NS2+MT)*EXP2 

307 CONTINUE 

c ET=ET+LEADTERM*(ET1+ET2) 

c EP=EP+LEADTERM*(EP1+EP2) 

305 CONTINUE 

C**» ******************************************** ' end mode summation 
ETT=LEADTERM* ( ET1+ET2 ) 

EPT=LEADTERM* ( EP 1+EP2 ) 

******************************************add in antenna contribution 
C FEED CONTRIBUTION IN THE FAR FIELD IS ETF , EPF 
C *********************************************** *****BESSJ 1 . 

C USE CLOSED FORM FEED EXPRESSION FROM SUBROUTINE ANTFF IF 
C IFF=0 ; OTHERWISE USE BRUTE FORCE EVALUATION FROM CIRCRTP 
C 

R0BS=1000 . 

IF(IFF.EQ.O) THEN 

CALL ANTFF(THR,PHR,THS,PHS, ARAD, ETF, EPF) 

C snte(6,*) ’ ET = ’ ,CABS(ET) , ’ ETF = ’ .CABS(ETF) 

ELSE 

RHB=R0BS*SIN(THR) 

ZHB=R0BS*C0S(THR) 

CALL CIRCRTP (CNPHI , XP , AP , ARAD , THS , PHS , 

& PHR , RHB , ZHB , CIRCR , CIRCT , CIRCP ) 

C REMOVE THE l/R DEPENDENCE BECAUSE EXP(-jkR)/R IS OMITTED IN 
C THE SCATTERED FIELDS ET AND EP 
ETF=CIRCT*ROBS 
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c 



EPF=CIRCP*ROBS 

orite(6 , * ) ’ ET = ’,CABS(ET),' ETF(ROBS) = ’ .CABS(ETF) 

ENDIF 

C*** ***************************************** **total E theta and E phi 
C write (6,*) ’THR( 1 ,jj,j,ii,i,NM, * ) = ’,180*THR/pi 

ET CO MB = ETT + ETF 
EPCOMB=EPT+EPF 

EMAG=SQRT(CABS(ETC0MB)**2+CABS(EPC0MB) **2) 

EMAGNONE=SQRT(CABS(ETF) **2+CABS (EPF) **2) 

IF(EMAG . GT . EMAX) THEN 
EMAX=EMAG 
EMAXNONE=EMAGNONE 

WRITE(6 , *) 1 EMAG( 1 , J , I , 1 ) = * , EMAG 
END IF 

C****************************************det max field intensity 
UMAG-EMAG**2 
UMAGN0NE=EMAGN0NE**2 
IF (UH AG. GT. UMAX) THEN 
UMAX=UMAG 

THETAMAX=180*THR/PI 

PHIMAX=180*PHR/PI 

C WRITE(6 , *) 1 UMAX = \UMAX 

C WRITE(6 , *) 1 THETAMAX = ’ ,THETAMAX 

C WRITE(6 , *) ’ PHIMAX = » ,PHIMAX 

END IF 

IF(UMAGNONE . GT . UMAXNONE)THEN 
UMAXNONE=UMAGNONE 
END IF 

C***********************multiply integration, by beam solid angle (STRA) 
SUM=SUM+UMAG*STRA 
SUMNONE=SUMNONE+UMAGNONE*STRA 
301 CONTINUE 

£*************************************************- end theta integration ~ 

300 CONTINUE 

C********** t *** ******************************************* ~ end phi integration * 
C****** t******************** sca i e each summation by (b-a)/2 for Gauss integration 
PRAD=T1*P1*SUM 
PRADN0NE=T1*P1*SUMN0NE 

WRITE(6 ,*) ’PRAD = ’ , PRAD 

WRITE(6 ,*) ’PRAD NO RADOME = * , PRADNONE 

C******** c************************************************************* ********** 

C End Integration over closed surface radiation intensity * beam solid angle 

£*********** ********************************************************************* 
C* *************** ** ************************************** *********** compute gain 
GAIN=4*PI*UMAX/PRAD 
GAINN0NE=4*PI*UHAXN0NE/PRADN0NE 
GDB=10.*alogl0(GAIN) 

GDBN0NE=10 . *aloglO(GAINNONE) 

WRITE (6,*) ’GAIN ’ .NDIVT, ’ = ‘.GAIN,’ IN DB = ’ , GDB 

WRITE(6 ,*) ’GAINNONE = ’ .GAINNONE, ’ IN DB = ’ .GDBNONE 

WRITE (7,*) ’ NDIVT (’ .NDIVT, > )=’ .NDIVT, ’ ; NDIVP(’ .NDIVP, > )=- .NDIVP 
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WRITE (7 , 2003 ) ’ GAIN( ’ , NDIVT , ’ )=’ .GAIN, ’ ;GAINDB( ’ , NDIVT , > )=’ , GDB 
WRITE (7 , 2003 ) ’ GAINNONEC ’ , NDIVT , ’ ) = ’.GAINNONE, 

& ’ ; GAINNONEDB( ’ , NDIVT , ’ ) = ’ .GDBNONE 

2003 FORMAT ( / (A,I3,A,F12.5,A,I3,A,F10.5) ) 

C500 continue 
STOP 

999 end 

C************************************************************£ND op MAIN PROGRAM 
0***** ******************************************* **** SUBROUTINE PLANE 
C REFERENCE: TECHNICAL REPORT TR-80-1 ; (pages 57-64) 

SUBROUTINE PLANE ( Ml , M2 , NP , NT , RH , ZH , XT , AT , THR , R) 

C 

C PLANE WAVE EXCITATION VECTOR IN THE FAR FIELD. FROM MAUTZ AND HARRINGTON. 

C MODIFIED TO DO ONLY ONE ANGLE AND FREQUENCY PER CALL. 

C ** NEW BESSEL FUNCTION ROUTINE FROM NUM. RECIPES ** 

C 

COMPLEX R( 1600) , U,U1 ,UA ,UB ,FA( 1500) ,FB(1500) ,F2A ,F2B ,F1A ,F1B 
COMPLEX U2,U3,U4,U5 

DIMENSION RH (400 ) ,ZH(400) ,XT(4) ,AT(4) ,R2(4) ,Z2(4) 

MP=NP-1 
MT=MP- 1 
N=MT+MP 
N2=2*N 
CC=COS (THR) 

SS=SIN(THR) 

U=(0. ,1.) 

Ul=3. 141593*U**M1 
M3=M1+ 1 
M4=M2+3 

IF(Ml.EQ.O) M3=2 

M5=M1 +2 

M6=M2+2 

DO 12 IP=1 , MP 

K2=IP 

I=IP+1 

DR= . 5* (RH(I)-RH(IP ) ) 

DZ= . 5*(ZH(I)-ZH(IP) ) 

D1=SQRT(DR*DR+DZ^DZ) 

Rl= . 25* (RH(I)+RH(IP) ) 

IF(R1 . EQ . 0 . ) Rl=l. 

Zl= . 5*(ZH(I)+ZH(IP) ) 

DR= . 5*DR 
D2=DR/R1 
DO 13 L=1 , NT 
R2(L)=R1+DR*XT(L) 

Z2 (L)=Z1+DZ*XT(L) 

13 CONTINUE 
D3=DR*CC 
D4=-DZ*SS 
D5=D1*CC 
DO 23 M=M3,M4 
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FA(M)=0 . 

FB(H)=0. 

23 CONTINUE 

DO 15 L=1 ,NT 
X=SS*R2(L)*2 . 

ARG=Z2(L)*CC 

UA=AT(L)*CMPLX(COS( ARG) , SIN(ARG) ) 

UB=XT(L)*UA 

C THIS LINE REPLACES HARRINGTON'S BLOCK 
DO 25 M=M3 , M4 
BES=BESSJ (M-2 , X) 

FA(M)=BES*UA+FA(H) 

FB(M)=BES*UB+FB(M) 

25 CONTINUE 
15 CONTINUE 

IF(Hl.NE.O) GO TO 26 
FA(1)=-FA(3) 

FB(1)=-FB(3) 

26 UA=U1 

DO 27 H=H5 ,M6 
M7=H-1 
M8=H+ 1 

F2A=UA*(FA(M8)+FA(H7) ) 

F2B=UA*(FB(M8)+FB(M7) ) 

UB=U*UA 

F1A=UB*(FA(M8)-FA(M7) ) 

F1B=UB*(FB(M8)-FB(M7) ) 

U4=D4*UA 

U2=D3*F1A+U4*FA(M) 

U3=D3*FlB+U4*FB(M) 

U4=DR*F2A 

U5=DR*F2B 

K1=K2-1 

K4=K1+N 

K5=K2+N 

R(K2+MT)=-D5*(F2A+D2*F2B) 

R(K5+HT)=D1*(F1A+D2*F1B) 

IF(IP.EQ. 1) GO TO 21 
R(K1)=R(K1)+U2-U3 
R(K4)=R(K4)+U4-U5 
IF(IP.Eq.MP) GO TO 22 

21 R(K2)=U2+U3 
R(K5)=U4+U5 

22 K2=K2+N2 
UA=UB 

27 CONTINUE 
12 CONTINUE 

RETURN 

END 

C*** ******************************** ******SUBROUTINE ANTFF 
SUBROUTINE ANTFF (THETA , PHI , THETAS , PHIS , ARAD , ETF , EPF ) 
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C COMPUTES CLOSED FORM ANTENNA PATTERN IN THE FAR FIELD. 

C ETF AND EPF ARE RETURNED. ANGLES ARE IN RADIANS. 

C PRESENT ENTRY IS FOR: 

C 

C UNIFORM ILLUMINATION SCANNED TO A DIRECTION (THETAS , PHIS) 

C 

COMPLEX ETF , EPF , JK , SCL , EB 
JK=(0. 0,6. 283185307) 

BK=6. 283185307 
PI=BK/2 . 

C SET RADIATION PATTERN TO ZERO IN THE REAR HEMISPHERE 
IF(COS(THETA ) . LT . 0 . ) THEN 
EB=(0. ,0. ) 

ELSE 

BB=SIN (THETA) *SIN( PHI) -SIN (THETAS) *SIN (PHIS) 

AA=SIN( THETA )*COS( PHI ) -SIN (THETAS ) *COS( PHIS) 
ARG=BK*ARAD*SQRT(AA**2+BB»*2) 

SCL=-(0. , 1 . )*2. *PI**2 
IF(ABS(ARG).LT.1.E-5)THEN 
EB=SCL/2 . 

ELSE 

EB=SCL*(BESSJ1(ARG)/ARG) 

ENDIF 

ENDIF 

C EB IS THE X COMPONENT: GET ETHETA (ETF) AND EPHI (EPF) 

C ASSUME AN ELEMENT FACTOR, ELHT 
ELMT=SQRT(ABS(COS(THETA) ) ) 

ETF=EB*COS (THETA') *COS (PHI ) *ELMT 
EPF=-EB»SIN(PHI ) *ELHT 
RETURN 
END 

£**«******««*****«**«******* *****«*************3ESSEL FUNCTIONS 
FUNCTION BESSJ(N.X) 

C RETURNS THE BESSEL FUNCTION B OF ORDER N (>1) AND REAL 
C ARGUMENT X. 

PARAMETER (IACC=40 , BIGN0=1 . E10 ,BIGNI=1 . E-10) 

IF(N.EQ.O) THEN 
C IF N=0 CALL BESS JO 

BESSJ=BESSJO(X) 

ELSE IF(N.EQ.l) THEN 
C IF N=1 CALL BESSJ1 

BESSJ=BESS J1 (X) 

ELSE 

C IF N>1 USE RECURSION 
BESSJ=0 . 

IF(X . NE . 0 . ) THEN 
T0X=2./X 

IF(X . GT. FLOAT(N) ) THEN 
BJM=BESSJO(X) 

BJ=BESSJ1 (X) 

DO 11 J=1,N-1 
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B JP= J*TOX*B J-B JM 
B JM=B J 
B J=BJP 

11 CONTINUE 
BESS J=B J 
ELSE 

M=2*( (N+INT(SQRT(FLOAT(IACC*N) ) ))/2) 

BESS J=0 . 

JSUM=0 . 

SUM=0 . 

B JP=0 . 

B J = 1 . 

DO 12 J=M ,1,-1 
BJM= J*TOX*B J-B JP 
B JP=B J 
B J=B JM 

IF(ABS(BJ) .GT.BIGNO) THEN 

B J=B J+BIGNI 

B JP=B JP*BIGNI 

BESS J=BESS J*BIGNI 

SUM=SUM*BIGNI 

ENDIF 

IF(JSUM.NE.O) SUM=SUM+BJ 
JSUM=1-JSUM 
IF(J.EQ.N) BESS J=B JP 

12 CONTINUE 
SUM=2 . *SUM-B J 
BESS J=BESS J/SUM 
ENDIF 

ENDIF 

ENDIF 

RETURN 

END 

FUNCTION BESSJO(X) 

C 

C BESSEL FUNCTION OF 0 ORDER, REAL ARGUMENT X 
C (SEE ' NUMERICAL RECIPES', P.172) 

C 

REAL+8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6, 

* SI , S2 , S3 , S4 , S5 , S6 

DATA PI ,P2,P3,P4,P5/l .DO,-. 109862827D-2 , . 27345 10407D-4 , 

* - . 2073370639D-5 , . 2093887211D-6/ 

DATA Q1 , Q2 , Q3 ,Q4 ,Q5/- . 1562499995D-1 , . 1430488765D-3 , 

* -.6911 147651D-5 , . 7621095161D-6 , - . 934945152D-7/ 

DATA R1 ,R2 , R3 ,R4 , R5 , R6/57568490574 . DO , - 13362590354 . DO , 

* 651619640 . 7D0 , -11214424 . 18D0 , 77392 . 33017D0 , -184 . 9052456D0/ 
DATA S 1 , S2 , S3 , S4 , S5 , S6/5756849041 1 . DO , 1029532985 . DO , 

* 9494680 . 718D0 , 59272 . 64853D0 , 267 . 8532712D0 , 1 . DO/ 

IF(X . EQ .0 . ) THEN 

BESS J0=1 . 

ELSE IF(ABS (X) . LT . 8 . ) THEN 
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Y=X**2 

BESSJ0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6) ) ) ) )/ 
* (S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) ) ) ) ) 

ELSE 

AX=ABS(X) 

Z=8 . / AX 
Y=Z**2 



XX=AX-. 785398164 

BESSJO=SqRT( .636619772/ AX)* (C0S(XX)*(P1+Y*(P2+Y*(P3+ 

* Y*(P4+Y*P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+ 

* Y*(Q4+Y*Q5) ) ) ) ) 

ENDIF 

RETURN 

END 

FUNCTION BESSJl(X) 

C 

C BESSEL FUNCTION B OF ORDER 1, REAL ARGUMENT X 
C (SEE ’NUMERICAL RECIPES’, P.173) 

C 

REAL*8 Y ,P1 , P2 ,P3 ,P4 , P5 , Q1 , Q2 , Q3, Q4 ,Q5 ,R1 ,R2 ,R3 ,R4,R5 ,R6 , 

* S 1 , S2 , S3 , S4 , S5 , S6 

DATA PI , P2 ,P3 ,P4, P5/1 . DO , . 183105D-2 , - . 3516396496D-4 , 

* . 2457520 174D-5 , - . 203370 19D-6/ 

DATA Q 1 , Q2 , Q3 , Q4 , Q5/ . 04687499995D0 , - . 2002690873D-3 , 

* . 8449199096D-5 99228987D-6 , . 105787412D-6/ 

DATA R1 , R2 , R3 ,R4 , R5 , R6/72362614232 . DO , -7895059235 . DO , 

* 242396853 . IDO , -2972611 . 439D0 , 15704 . 4826D0 , -30 . 16036606D0/ 
DATA SI ,S2 , S3 ,S4 , S5 , S6/144725228442 . DO , 2300535178 . DO, 

* 18583304 . 74D0 , 99447 . 43394D0 , 376 . 9991397D0 , 1 . DO/ 
IF(X.Eq.O.) THEN 

BESS J1=0 . 

ELSE IF(ABS(X) . LT . 8 . ) THEN 
Y=X**2 

BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6) ) ) ) )/ 

* (Sl+Y* (S2+Y*(S3+Y*(S4+Y*(S5+Y*S6) ) ) ) ) 

ELSE 

AX=ABS(X) 

Z=8./AX 
Y =Z**2 



XX=AX-2. 356194491 

BESSJ1=SQRT( .636619772/AX)* (COS (XX )*(P1+Y*(P2+Y*(P3+ 

* Y*(P4+Y*P5) ) ) )-Z*SIN(XX)*(qi+Y*(q2+Y*(Q3+ 

* Y*(q4+Y*q5)))))*SIGN(l . ,X) 

ENDIF 

RETURN 

END 
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c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c***************************************************************************** 



COMMENTS : THIS SUBROUTINE HILL RIGOROUSLY CALCULATE THE ELECTRIC FIELD 

FOR A CIRCULAR APERTURE. THE FIELD IS CALCULATED AT THE COORDINATES 
SPECIFIED BY RH( . ) AND ZH(.) THE APERTURE IS LOCATED AT Z = 0 , 

AND SCANNED TO A DIRECTION (THS.PHS). THE SUBROUTINE OGIVE IS THE 
SOURCE OF THE GEOMETRIC DATA REQUIRED BY CIRC TO PERFORM COMPUTATIONS. 
ALL PHYSICAL DIMENSIONS ARE NORMALIZED TO WAVELENGTH. 
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c***************************************************************************** 
SUBROUTINE CIRCRTP( CNPHI , XP , AP , ARAD , THS , PHS , 

* PHI ,RZ ,Z ,CIRCR, CIRCT ,CIRCP) 

INTEGER CNPHI 
DIMENSION XP(IOO) ,AP(100) 

REAL K 

COMPLEX J, JK,G1 , G2,X1 ,CON , CC,SUMX,SUMY, SUMZ 
COMPLEX CIRCR, CIRCT, CIRCP 
PI=3 . 141592654 
ETA=120 . *PI 

C SET FIELD COMPONENTS TO ZERO IN THE REAR HEMISPHERE 
CIRCR=(0 . ,0. ) 

CIRCT=(0. ,0. ) 

CIRCP=(0 . ,0. ) 

C Z OF THE OBSERVATION POINT IS >= 0 IN THE FORWARD HEMISPHERE 
IF(Z.GE.O. ) THEN 
SUMX=(0. 0,0.0) 

SUMY=(0. 0,0.0) 

SUMZ=(0. 0,0.0) 

K=6. 283185307 
J = (0.0, 1.0) 

JK= (0.0,6.283185307) 

C0N=-J/(4.*PI) 

C OMIT THE FACTOR OF ETA SINCE IT IS NOT INCLUDED IN SUBROUTINE ZMAT . 

C ARAD**2 IS FOR SCALING OF GAUSSIAN INTEGRATION. READJUST SCALING 

C TO AGREE WITH OLD SUBROUTINE SPHERE (ARBITRARY) 

CC= CON 

C OUTER LOOP: INTEGRATE OVER X, -ARAD < X < ARAD. 

C INNER LOOP: INTEGRATE OVER Y, -ARAD < Y < ARAD. 

C EVALUATE ANTENNA FIELD AT RZ,Z: 

XO=RZ*COS(PHI) 

YO=RZ*SIN(PHI) 
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STS=SIN (-THS ) 

S=RZ/SQRT(RZ**2+Z**2) 

C=Z/SQRT(RZ**2+Z**2) 

C INTEGRATE IN X,Y COORDINATES 
DO SO M=1 ,CNPHI 
X=ARAD*XP(M) 

DO 60 N=1,CNPHI 
Y=ARAD*XP(N) 

RP=SQRT(X**2+Y**2) 

C SKIP INTEGRATION POINTS BEYOND THE ANTENNA RADIUS 
IF(RP.LE. ARAD) THEN 
PHIPRIME=ATAN2 (Y , X+l . E- 10) 

R=SQRT( ( (RZ-RP ) **2+Z**2)+(4*RZ*RP*SIN( (PHI-PHIPRIME)/2)**2) ) 
G1=(((K*R)**2)-1-(JK*R))/R**3 
G2=(3+(3* JK*R)-(K*R) **2)/R**5 
X1=JK*(RP*C0S(PHIPRIME-PHS)*STS-R) 

Y1=X0-X 

X2=Y1**2 

Y2=Y0-Y 

SUMX=SUMX+(CC*AP(M)*AP(N)*(G1+(X2*G2))*CEXP(X1) )*TAPER(RP) 
SUMY=SUHY+(CC*AP(M)*AP(N)*YH'Y2*G2*CEXP(X1))*TAPER(RP) 
SUMZ=SUMZ+(CC*AP(M) *AP(N) *Y1*Z*G2*CEXP(X1 ) ) *TAPER(RP) 

ENDIF 

60 CONTINUE 

SO CONTINUE 

C ASSUME AN ELEMENT FACTOR, ELMT 
ELMT=SQRT(ABS(C) ) 

CIRCT=(C*COS(PHI)*SUMX+C*SIN(PHI)*SUMY-S*SUMZ)*ELMT 

CIRCR=(S*COS(PHI)*SUMX+S*SIN (PHI)*SUMY+C*SUMZ)*ELMT 

CIRCP=(-SIN(PHI)*SUMX+COS(PHI)*SUMY)*ELMT 

ENDIF 

RETURN 

END 

c *********************************** *******function taper 

FUNCTION TAPER(RHO) 

C SPECIFY ANTENNA ILLUMINATION FUNCTION. REAL FUNCTION OF 
C RHO ONLY 

TAPER=1 . 

RETURN 

END 
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